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Abstract 

We review recent progress toward establishing lattice Quantum Chromodynamics as a predictive 
calculational framework for nuclear physics. A survey of the current techniques that are used to 
extract low-energy hadronic scattering amplitudes and interactions is followed by a review of recent 
two-body and few-body calculations by the NPLQCD collaboration and others. An outline of the 
nuclear physics that is expected to be accomplished with Lattice QCD in the next decade, along 
with estimates of the required computational resources, is presented. 

1 Introduction 

It has been known for the past four decades that Quantum Chromodynamics (QCD), together with the 
elect roweak interactions, underlies all of nuclear physics, from the s-wave nucleon-nucleon scattering 
lengths to the mechanism of synthesis of heavy elements in stars. Soon after the discovery of QCD it 
became apparent that the complexities of the theory at strong coupling would hinder analytic progress 
in understanding the properties of the simplest hadrons, let alone the simplest features of the nuclear 
forces. Wilson pointed the way to eventual direct quantitative confirmation of the origins of nuclear 
physics by formulating Lattice QCD P, a regularization and non-perturbative definition of QCD that 
is suitable for the intensive computational demands of solving QCD in the infrared. 

Only in the last five years or so has Lattice QCD emerged from a long period of research and develop- 
ment — where only qualitative agreement between calculations and experiments could be claimed — to 
the present, where precise predictions for hadronic quantities are being made. In particular, presently, 
fully-dynamical calculations with near-exact chiral symmetry at finite lattice-spacing have become stan- 
dard, with lattice volumes of spatial extent 2.5 fm and with lattice spacings in the range 0.12 fm. 
Very recently, preliminary calculations at the physical quark masses have been carried out. However, 
it is still the norm that the light-quark masses, m^, are larger than those of nature, with typical pion 
masses 171^^ ^ 300 MeV. It is expected that in the next five years, calculations at the physical light- 
quark masses, ^ 140 MeV, in large volumes, 6 fm, and at small lattice spacings, 0.06 fm 
will become commonplace. 

Nuclear physics is a vast, rich field, whose phenomenology has been explored for decades through 
intense experimental and theoretical effort. However, there is still little understanding of the connection 
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to QCD and the basic building blocks of nature, quarks and gluons. For instance, as a first "benchmark- 
ing" step. Lattice QCD should post-diet the large nucleon-nucleon scattering lengths and the existence 
of the deuteron, the simplest nucleus. The connection between QCD and nuclear physics will be firmly 
established with Lattice QCD, and will allow for an exploration of how nuclei and nuclear interactions 
depend upon the fundamental parameters of nature. In particular, it is believed that an understanding 
of the fine-tunings that permeate nuclear physics will finally be translated into fine-tunings of the light- 
quark masses. While these issues are of great interest, and it is important to recover what is known 
experimentally to high precision, these goals are not the main objective of the Lattice QCD effort in 
nuclear physics. The primary reason for investing resources in this area is to be able to calculate phys- 
ical quantities of importance that cannot be accessed experimentally, or which can be measured with 
only limited precision in the laboratory. Two important examples of how Lattice QCD calculations can 
impact nuclear physics are in the structure of nuclei and in the behavior of hadronic matter at densities 
beyond that of nuclear matter. 

Refined many-body techniques for studying the structure of nuclei, such as Greens function Monte- 
Carlo (GFMC) [2j, already exist. This method has led to calculations of the ground states and excited 
states of light nuclei, with atomic number 14. Using only the modern nucleon-nucleon (NN) poten- 
tials that reproduce all scattering data below inelastic thresholds with x^/dof ^ 1, such as AVis f3], 
these models fail to recover the observed structure of light nuclei. However, the inclusion of a three- 
nucleon interaction greatly improves the predicted structure of nuclei [2j. Lattice QCD will be able 
to calculate the interactions of multiple nucleons, bound or unbound, in the same way that it can be 
used to determine the two-body scattering parameters. For instance, a calculation of the three-neutron 
interactions will be possible. 

One of the great challenges facing nuclear physics is to determine the behavior of hadronic matter at 
densities away from that of nuclear matter, such as those that occur in the interior of neutron stars [4]. 
There are a number of possibilities for the composition of the hadronic matter at the center of a neutron 
star. One possibility is that it is composed of neutrons and protons. Another possibility is that it is 
composed of neutrons, protons and a kaon condensate, and a third possibility is that it is composed of 
neutrons, protons and S~'s. At present, all of these compositions (and other more exotic scenarios) are 
possible theoretically, and not excluded experimentally or observationally. On the theoretical side, the 
main uncertainty in establishing the composition of hadronic phases, and hence the equation of state, is 
the interactions between the strange hadrons, such as the kaon and S~, and the protons and neutrons. 
These interactions are poorly known experimentally due to the short lifetime of the strange hadrons, 
however Lattice QCD promises to reliably calculate these interactions from QCD with quantifiable 
uncertainties. Indeed lattice QCD calculations will provide the best determinations of hyperon-nucleon 
and kaon-nucleon scattering amplitudes and hence will play a crucial role in determining the role of 
hyperons and strange mesons in neutron stars. 

Lattice QCD is a technique in which Euclidean space correlation functions are calculated by a 
Monte-Carlo evaluation of the Euclidean space path integral [Ij. The calculations are performed in 
Euclidean space so that field configurations that have a large action are exponentially suppressed. This 
is in contrast with Minkowski space in which large action contributions result in a complex phase 
which will average to an exponentially small contribution with nearby configurations. In this approach, 
space-time is discretized with the quarks residing on the lattice sites, and the gluon fields residing on 
the links between lattice sites. The lattice spacing, 6, the distance between adjacent lattice sites, is 
required to be much smaller than the characteristic hadronic length scale of the system under study. The 
effects of a finite lattice spacing can be systematically removed by combining calculations of correlation 
functions at several lattice spacings with the low-energy effective field theory (EFT) which explicitly 
includes the discretization effects. This type of EFT is somewhat more complicated than its continuum 
counterpart as it must reproduce matrix elements of the Symanzik action constructed with higher 
dimension operators induced by the discretization [5]. While the action lacks Lorentz invariance and 
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rotational symmetry, it is constrained by hypercubic symmetry. As computers have finite memory 
and performance, the lattice volumes are finite in all four space-time directions. Generally, periodic 
boundary conditions (BC's) are imposed on the fields in the space-directions (a three-dimensional torus), 
while (anti-)periodic BC's are imposed on the (quark) gauge fields in the time-direction, which in many 
cases is much larger than the space-directions (in order to approach the zero-temperature limit). 

For the calculations we will be discussing in this review, the lattice volumes are large compared 
with the Compton wavelength of the pion, and deviations of single-particle properties from their infinite 
volume values are exponentially small, generically ^ ^-m^L ^ These finite- volume eflFects may be removed 
by finite- volume EFT. Perhaps the most important "economic" feature of Lattice QCD calculations is 
that the computational resources required to perform lattice calculations increase with decreasing quark 
mass, and presently, no calculations exist at the physical quark masses, ttIt^ ^ 140 MeV, that have 
m^T-L 3> 1. It remains the case that lattice calculations are performed at unphysical values of the quark 
masses, and the light quark mass dependence of the observable of interest, which can be determined 
perturbatively in the low-energy EFTs, is used to extrapolate to the physical light quark masses. 
Therefore, the practical situation with current Lattice QCD calculations is that they are performed at 
finite lattice spacing, within finite volumes and at unphysical quark masses. The appropriate EFT (e.g. 
chiral perturbation theory (x-PT), heavy-baryon-^-PT (HBx-PT) ) is then used to extrapolate to the 
infinite volume, continuum limit of QCD where physical predictions can be made. 



2 Lattice QCD Technology 

Quantum Chromodynamics (QCD) can be defined non-perturbatively as the continuum limit of a 
Lattice gauge theory. This approach provides both an ultraviolet regulator of the continuum field 
theory and admits numerical evaluation of the functional integrals required for calculating physical 
observables. In the continuum, the QCD path integral is[^ 

Z = J VA^V^V^ e-^'^''' ^^"^'^ 
/ 

where /^qcd is the QCD Lagrange density, is the QCD gauge field (describing the gluons), is 
the gauge field strength, m/ is the quark mass, and ^0/, '0/ are the fermion fields representing the quark 
flavors. Dy, is the covariant derivative which ensures gauge invariance and 7^ are the Dirac-matrices. 
Explicitly, the covariant derivative acting on a quark field is given by 

Dy^{x) = dy^l;{x)+igsAy{x)^l;{x) , Ay{x) = A;{x) , (2) 

where = with the being the Cell-Mann matrices. The strong coupling gg appearing in 

eq. ([2]) is related to the strong fine-structure constant, a^, via ag = 5'^/(47r). The field strength tensor 
is defined in terms of the gluon field through 

Gy,{x) = dyA,{x) - d,Ay{x) + ig^ [ Ay{x) , A,{x) ] , G^,(x) = G^,(x) , (3) 

Observables (physical quantities) in this theory can be calculated from correlation functions of operators 
O that are functions of the quantum fields (quarks and gluons). 

{0) = ^ J VAyV^V^ O e-^^'^ ^Q^^ . (4) 



very nice introduction to Lattice QCD can be found in Ref. 
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Figure 1: A two dimensional slice of the four dimensional space-time lattice, /x and v denote 
unit-vectors in the indicated directions. ip{x) denotes a fermion-field at the lattice-site x, 
U^{x) denotes the gauge link from the lattice-site x to the site x + and P^y{x) denotes 
the 1x1 Wilson plaquette centered at x + 6/x/2 + 6z//2. 



The path integral can be rigorously defined on a discrete space-time. In order to preserve gauge 
invariance the gauge fields are discretized as special unitary matrices, in the group aS[/(3), on the links 
of the lattice (see Figure [T]). The discrete gauge action is the sum over all plaquettes P^zy(x) which 

are the product of the links = exp J^^^ rfx'A^(x')^ around the elementary plaquettes of the 

lattice, 

^.(f^) = /9E(l-^^^TrP^,(x)') , (5) 

with 

Pi^u = U^{x)U,{x + ji)Ul{x + v)Ul{x) . (6) 

/3 is the lattice gauge coupling that is related to the strong coupling via /3 = 2Nc/ gl where Nc is the 
number of colors. Taking the naive continuum limit, this action becomes the familiar continuum gauge 
action, — J d^x\ {G^^^{x)) . The action in Eq. 5 is the well known Wilson gauge [Ij action, and while 
this discretization is not unique, it is the simplest. It can be modified by adding larger loops with 
coefficients appropriately chosen to achieve better convergence to the continuum, which is the ultimate 
goal of the calculation. 
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Inclusion of the quarks (fermions), which are defined on the vertices of the lattice, is a challenging 
problem. A naive discretization of a single continuum fermion field introduces sixteen light lattice 
fermion flavors in the four dimensions. In QCD only three light quarks (up, down and strange) and 
three heavy quarks (top, bottom and charm) are needed. The doublers, i.e. the additional light 
lattice fermion degrees of freedom, can be avoided by several ingenious formulations of lattice fermions. 
Wilson fermions, which were introduced first [Ij, eliminate the doublers by adding irrelevant dimension 
five operators in the action that lift the masses of the doublers, leaving only one light fermion in the 
spectrum. However, this explicitly breaks chiral symmetry and introduces lattice artifacts that scale 
as 0{b). The Wilson action can be improved by the addition of one dimension-5 bilinear operator, 
Osw = i^^^^Gj^jyip^ the Sheikholeslami-Wohlert term [7J, with a coefl&cient, Csw^ that can be tuned so 
that the lattice artifacts are parametrically reduced to scale as 0(6^). In addition, the lattice spacing 
artifacts can be reduced to 0(6^) by the addition of a twisted-mass term to the action [8j. The isospin 
violation that the twisted-mass term introduces into the action can be removed using chiral perturbation 
theory. Kogut-Susskind fermions [9j (staggered fermions) provide another way to remove some of the 
doublers and re-interpret the remaining four as four degenerate fiavors. In this approach a [/(I) chiral 
symmetry still remains unbroken and lattice artifacts scale as 0{P). Kogut-Susskind fermions become 
problematic when the required number of fiavors is not a multiple of four (as is the case for QCD). 
In addition, the broken fiavor and chiral symmetries introduce large lattice artifacts, although they 
scale as 0(6^). Finally, the so called domain wall fermions pUl [TTl [T2] and overlap fermions p!3l [T4] 
are fermionic actions that both preserve a lattice chiral symmetry at finite lattice spacing (they satisfy 
the Ginsparg- Wilson relation [15j 75!^ + D^^ = b D^^D) and are doubler free. Unfortunately, such 
formulations are computationally significantly more expensive. It should also be stated that other 
actions are being explored that approximately satisfy the Ginsparg- Wilson relation, for instance the 
fixed-point action [16], or the chirally-improved action, e.g. Ref. pTlfTS] . In all cases the lattice fermion 
action is of the form 

Sf = i,D{U)ij (7) 

where ifj is the fermion "vector" and D{U) is a sparse matrix acting on the fermion vector, that 
depends on the gauge field U. As an example of the form of Sf for a lattice action, its explicit form for 
the naive Wilson action is 

^Wilson = H^hAM^)H^ + bl^) - Ul{x-bfi)^{x-bfi)] + mfJ2 H^m^) ■ (8) 

X X 

The partition function in the case of two quark fiavors is 

/2,X X 

= jl[dU,{x) det{D{UyD{U)) e-'^^''^ (9) 



fl^X 



The integration over the quark fields, which are represented by Grassmann numbers, can be done exactly. 
In addition, the quark matrix D{U) represents one fiavor, however since detD{U)'^ = detD{U)^ the 
determinant det (^D{Uy D{U)) represents two fiavors. In the case of correlation functions, integrating 
out the quarks gives the following expression 

= 1/ UdU.i^^ ^^D{uy^^ {D{UYD{U)) e-'^^^^ , (10) 



^In certain cases, such as with overlap fermions, the matrix is not sparse but has sparse like properties, i.e. the matrix 
vector multiplication is a computationally "cheap" operation. 
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where the quantities O depend on the inverse of the quark matrix and possibly exphcitly on the gauge 
fields. The above expressions are only valid in the case of two fiavors of quarks (the up and the down) 
which both have the same mass, which is a good approximation to the low energy physics of QCD. A 

1 /9 

strange quark can be easily added by including det (^D{Uy D{U)) in the partition function. 



The computation of Eq. (10) is the main numerical task faced in Lattice QCD calculations. The 



integral over the gauge fields in Eq. (10) is of extremely large dimensionality. Given that QCD has a 
fundamental length scale of ^ 1 fm (10~^^ cm), calculations must be performed in lattice volumes that 
have a physical size much larger than 1 fm in order to control finite volume eflFects, and with lattice 
spacings much smaller than 1 fm in order to be close to the continuum limit. With moderate choices for 
the volume and the lattice spacing, a lattice volume of ^ 32^ x 256 is currently practical. With the color 
and spin degrees of freedom, such calculations involve ^ 10^ degrees of freedom. The only practical 
way that this type of computation can be done is by using Monte-Carlo integration. Fortunately, the 
combination of the quark determinant and the gauge action, 

V{U) = \ det {D{IJ)'^D{IJ)) e-^^(^) , (11) 

is a positive definite quantity which can be interpreted as a probability and hence importance sampling 
can be employed. The basic algorithm is to produce gauge field configurations {[/} with probability 
distribution V(JJ) and then evaluate 

1 ^ 1 



At finite A, the estimate of O is approximate, with an uncertainty that converges to zero as (9(1/ v A) 



Both for the gauge field configuration generation and the evaluation of Eq. (12), the linear system of 
equations 

Dt([/)[m]D([/)[m]x = </>, (13) 

needs to be solved where the dependence of the quark matrix on the quark mass m is made explicit. 
Since the quark matrix is sparse, iterative solvers such as conjugate gradient (CG) can be used. The 
condition number of the quark matrix is inversely proportional to the quark mass. Since the physical 
quark masses for the up and down quarks are quite small, the quark matrix has a large condition 
number. With current computer resources this linear system cannot be solved exactly at the physical 
quark mass. For that reason the calculation is performed at heavier quark masses and then extrapolated 
to the physical point. The vast majority of the computer time used in these calculations is devoted to 
the solution of this linear system both in the context of gauge field generation and in the later stage of 



the calculation of physical observables through Eq. (12). 

Realistic lattice calculations require quark masses that result in pion masses below 400 MeV, allowing 
chiral EFT's to be used with some reliability. In addition, a dynamical strange quark is required in 
order to guarantee that the low energy constants of the EFT determined with lattice QCD match those 
of the physical theory. Although this task seems formidable, in the last several years there have been 
developments that make phenomenologically interesting calculations now possible. 

The works that we will be reviewing involves the computation of correlation functions with fixed 
particle number, for instance, correlation functions with the quantum numbers of the deuteron, or the 
correlation functions with the quantum numbers of twelve charged kaons. As the number of particles 
increases, the complexity of forming the correlation function grows dramatically, due to the need for 
the quarks to satisfy the Pauli principle. An alternate way to describe systems with large numbers of 
baryons, for instance, is to introduce a baryon-number chemical potential into the partition function in 
Eq. ([T]). However, the evaluation of the partition function with such a chemical potential is complicated 
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by the fact that the measure of the integration is no longer positive definite, that is to say, it suffers from 
the weU-known sign-problem. Incremental progress has been made in dealing with the sign-problem, 
and a conceptual break through will be required in order to calculate the properties of systems at high 
densities and low-temperatures. Naively it appears that there should be a direct connection between the 
sign-problem and the behavior of multi-baryon correlation functions. However, to date, little connection 
between the two approaches to the numerical evaluation of the properties of finite-density systems has 
been established. A recent review of the status of calculations at finite density using chemical potentials 
can be found in Ref. [19]. 

2.1 Lattice Fermion Actions 

The emergence of fermions that respect chiral symmetry [TOl [TTl [T2l [T3l [T4] on the lattice was one of 
the major recent developments in Lattice QCD. These formulations of lattice fermions allow for the 
reduction of the lattice spacing errors and approach the continuum limit in a smooth manner. However, 
the computational resources required to calculate with these fermions are an order of magnitude larger 
than any other variant of lattice fermions. In addition, the development of improved Kogut-Susskind 
fermion actions [^O'^'^T] that significantly reduce the 0(6^) errors, allowed for "cheap" inclusion of quark 
loop eflFects in the QCD correlation functions computed on the lattice. With this formulation, volumes 
with spatial extent as large as L ^ 5.8 fm are being used with light-quark masses as low as l/20th 
of the strange quark mass depending on available computing resources. However the fact that Kogut- 
Susskind fermions represent four fiavors of quarks complicates calculations when two or one fiavors are 
needed. From the operational point of view the problem is solved by introducing the Kogut-Susskind 
determinant raised to the n//4 power (rooted) into the path integral, where n/ is the desired number of 
fiavors. The non-integer power of the quark determinant introduces non-localities in the lattice action. 
It has been argued that the long distance physics that survives the continuum limit is not affected by 
such non-localities [22l [23l Ell [25l [26l [27^ I^i addition, at finite lattice spacing, the pathologies arising in 
the Kogut-Susskind fermion formulation can be dealt with in staggered x-PT [23l EU [28] . Although no 
rigorous proof exists, empirical evidence indicates that Kogut-Susskind fermions do describe the correct 
physics as long as the continuum limit is taken before the chiral limit [27j The results presented in 
this review that are based upon mixed-action calculations on the MILC lattice ensembles assumes that 
the continuum limit of such calculations is QCD. 

Finally, significant simulation algorithm developments have been achieved that have sped up the 
gauge field generation by orders of magnitude. Developments in the integrators needed for the molecular 
dynamics [30j that form the core of the Hybrid Monte-Carlo algorithm |^ are one important com- 
ponent of the improvements. A second important component of the new algorithms is preconditioned 
Hybrid Monte-Carlo p3ll34ll35j and rational hybrid Monte-Carlo (RHMC) [38j for calculations involving 
an odd number of fiavors. Together with multiple time scales ^36| i37j for evolving different parts of the 
molecular dynamics Hamiltonian, calculations close to or at the physical quark masses became possible 
for a variety of lattice fermion actions including Wilson, Kogut-Susskind and domain- wall. For a status 
review of the current dynamical simulation algorithms the reader is referred to [39] . The result of these 
algorithmic developments was that calculations with Wilson or improved Wilson fermions again became 
feasible. In addition dynamical simulations with domain wall fermions are also possible at relatively 
small additional "cost" compared to "cheap" fermion variants at the same physical parameters. The 
choice of the computational approach is made based on the resources each collaboration has available 
as well as the physics objectives. 

^It should be noted that there are some members of the lattice community who believe that the rooted-staggered 
action is fundamentally flawed and its continuum limit does not correspond to QCD (for a summary of these arguments, 
see Ref. [W). We disagree with these arguments, however we acknowledge that there is no proof that the continuum 
limit of the rooted-Kogut-Susskind action corresponds to QCD. 



7 



2.1.1 Mixed Action Calculations 



The mixed-action calculations discussed in this review employed Kogut-Susskind fermions to represent 
the QCD vacuum polarization effects associated with the two light flavors (up/down quarks) and the 
somewhat heavier strange quark. This is done by using gauge conflgurations generated with the appro- 
priate Kogut-Susskind fermion determinants incorporated into the probability distribution that enters 
the path integral. Since this part of the computation is separated from the calculation of correlation 
functions, gauge flelds generated by other collaborations can be used. NPLQCD made use of a number 
of ensembles of gauge conflgurations generated by the MILC collaboration [40]. Domain wall fermions 
were used to describe all external (valence) quarks. Because of the chiral symmetry that domain wall 
fermions satisfy, all correlation functions satisfy chiral Ward identities, ensuring that the leading or- 
der (LO) chiral behavior is continuum-like. The small corrections appearing due to Kogut-Susskind 
fermions in the vacuum loops can be taken care of systematically in x-PT [HI [42l [43]. Compared 
to calculations with Kogut-Susskind fermions in the valence sector, this formulation results in better 
control of the chiral behavior and possibly smaller discretization errors and also simplifles calculations 
with baryons. This approach was flrst introduced by the LHP collaboration for the study of nucleon 
structure [44, |45l [46l |47l |48]. Because the valence and sea quark actions are different, such calculations 
are inherently partially quenched and are not unitary. This type of calculation is sometimes referred 
to as the mixed action scheme. Unlike conventional partially quenched calculations, which become 
unitary when the valence quark mass is tuned to the sea quark mass, unitarity cannot be restored 
by tuning the valence quark mass. The next best option is to tune the valence quark mass in such a 
way that the resulting pions have the same mass as those made of the sea Kogut-Susskind fermions. 
In this case unitarity should be restored in the continuum limit, where the Uf = 2 + 1 staggered ac- 
tion has an SU (12) l ® SU (12) ji U{l)v chiral symmetry due to the four-fold taste degeneracy of 
each flavor, and each tt, K or has 15 degenerate additional partners. At flnite lattice spacing this 
symmetry is broken and the taste multiplets are no longer degenerate, but have splittings that are 
0{a'^b'^) |20l EK HHl Eni [5l]. The domain wall fermion mass is tuned to give valence pions (kaons) that 
match the Goldstone Kogut-Susskind pion (kaon) This choice gives pions that are as light as possible, 
resulting in better convergence of the x-PT extrapolation of the lattice results to the physical quark 
mass point. This tuning was used by the LHPC collaboration |44[ |45l [46l |47l [521 [53] . 



2.1.2 Anisotropic Clover Wilson Fermions 

Recently anisotropic lattices have proven useful for spectroscopy projects, and as the calculations needed 
for studying multi-hadron systems rely on spectroscopy, NPLQCD has recently adopted clover-improved 
Wilson fermion actions. In particular, the rif = 2 + 1 flavor anisotropic Clover Wilson action [541 [55] 
with stout-link smearing jS] of the spatial gauge flelds in the fermion action with a smearing weight 
of p = 0.14 has been used. The gauge flelds entering the fermion action are not smeared in the 
time direction, thus preserving the ultra-locality of the action in the time direction. Further, a tree- 
level tadpole-improved Symanzik gauge action with no 1 x 2 rectangle in the time direction is used. 
Anisotropy allows for a better extraction of the excited states as well as additional confldence that 



plateaus in the eflFective mass plots formed from the correlation functions (discussed in Section 2.2) 
have been observed, signiflcantly reducing the systematic errors in observables due to fltting. The 
gauge fleld generation is done by the Hadron Spectrum Collaboration (HSC), and these gauge fleld 
conflgurations have been used for excited hadron spectrum calculations by HSC [57i 



^This is the only Goldstone boson that becomes massless in the chiral limit at finite lattice spacing. 
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2.2 Euclidean Space Correlation Functions 



Most Euclidean space correlation functions computed in LQCD calculations (suitably Fourier trans- 
formed) are the sums of exponential functions. The arguments of the exponentials are the product of 
Euclidean time with the eigenvalues of the Hamiltonian associated with eigenstates in the finite- volume 
that couple to the hadronic sources and sinks. For a lattice that has infinite extent in the time-direction, 
the correlation function at large times becomes a single exponential dictated by the ground state en- 
ergy and the overlap of the source and sink with the ground state. As an example, consider the pion 
two-point function, C^+{t)^ generated by a source (and sink) of the form 7r+(x, t) = t^(x, t)75rf(x, t), 

CMt) = J](0|7r-(x,t)7r+(O,0) |0) , (14) 

X 

where the sum over all lattice sites at each time-slice, t, projects onto the p = spatial momentum 
states. The source 7r+(x, i) not only produces single pion states, but also all states with the quantum 
numbers of the pion. More generally, the source and sink are smeared over lattice sites in the vicinity 
of (x, t) to increase the overlap onto the ground state and lowest-lying excited states. Translating the 
sink operator in time via 7r+(x, t) = e^V+(x, 0)e~^^ and inserting a complete set of states, gives |^ 

= 5](0|7r-(x,0)|n)(n|7r+(O,0)|0) ^ • (15) 

At finite lattice spacing, the correlation functions for Wilson fermions remain sums of exponential 
functions, but for particular choices of parameters used in the domain- wall discretization, the correlation 
functions exhibit additional sinusoidally modulated exponential behavior at short-times with a period 
set by the lattice spacing [i6JJ . 

It is straightforward to show that the lowest energy eigenvalue extracted from the correlation function 



in Eqs. (14) and (15) corresponds to the mass of the 7r+ (and, more generally, the mass of the lightest 
hadronic state that couples to the source and sink) in the finite volume. The masses of stable single 
particle states can be extracted from a Lattice QCD calculation with high accuracy as long as the lattice 
spatial extent is large compared to the pion Compton- wavelength |^ 

Once a correlation function is calculated, a common objective is to extract the argument of the 
exponential function that persists at large times. One way to do this is to simply fit the function over a 
finite number of time-slices to a single exponential function. A second method, that is somewhat more 
useful in visually assessing the quality of the calculation, is to form the eflFective mass (EM) function, 
e.g 

M^^it'^tj) = 1 logf-%^%-) ^ , (16) 



where both t and Meff.(t;tj) are in lattice units. At large times, MQ^,{t]tj) becomes a constant equal 
to the mass of the lightest state contributing to the correlation function The anti-periodic boundary- 
conditions in the time-direction, imposed on the quark-fields in order to recover the correct fermionic 
partition function, result in the single meson correlation functions being periodic in the time direction 



^We assume the absence of external electroweak fields that may exert forces on hadrons in the lattice volume. 
^Finite- volume effects are exponentially suppressed [62 by factors of e~^^^. 

''This is obviously the most simphstic approach to this problem. One well-known method to extract the ground state 
and excited state energies is that of Liischer and Wolff [63l [65] in which the correlation functions resulting from different 
sources and sinks are calculated. The resulting matrix of correlation functions is diagonalized, and the EM function for 
each resulting eigenvalue can be used to extract the spectrum. 
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and symmetric about the mid-point in the time direction. They are the sum of cosh- functions, and for 
such correlation functions it is useful to construct the EM using 



Meflp.(t;tj) = — cosh ( ^ ^ ^^^^ 1 , (17) 

which becomes a constant value for a single exponential or a single hyperbolic-cosine (near the mid-point 
of the time-direction of the lattice). 



2.3 Statistical Analysis Methods 

Since Monte-Carlo integration is used to compute the relevant correlation functions, the statistical un- 
certainty must be carefully determined. The main observables extracted from the calculations presented 
in this review are energy eigenvalues and their differences, which contain information about phase shifts, 
scattering lengths and the three body interaction. The extraction of energy eigenvalues is done by fit- 
ting the relevant correlation functions to a sum of exponentials (or hyperbolic cosine functions). The 
optimal values for the energy are extracted from correlated x^-minimization fits that take into account 
the time correlations in the lattice calculations. In particular, the relevant parameters, such as the 
energies and the amplitude of each state that contributes to the correlation function, are determined as 
those that minimize 

X\A) = J2 [G{U) - A)] C-/ [G{tj) - F{tj, A)] , (18) 

where G{t) are the lattice two point correlation functions, F(t, ^4) are the fitting functions used, A 
denotes the set of fitting parameters over which x^{A) is minimized, and Cij is the covariance matrix. 
The lattice two point correlation functions are determined as averages over Monte-Carlo samples of 
the correlation function, Gk{t): 

1 

G{t) = -YG.it) , (19) 



N ^ 

k=i 



and 

N 



(20) 



The (standard) errors on the fitted parameters are determined by the boundaries of the error ellipsoid]^ 
In computing scattering parameters, the procedure for determining the statistical uncertainties is 
somewhat more involved due to the highly non-linear relation between the scattering amplitude and 
the energy levels of the two-hadron system. First, one is interested in the energy diflFerences between 
the energy levels of the two-hadron system and the sum of the masses of the two free hadrons (similarly 
for the case of three or more hadrons). These energy differences can be determined in two ways. 
The simplest is where ratios of correlators are constructed in such a way so that they are a sum over 
exponentials parametrized by the desired energy splittings. In this case Jackknife is used to determine 
the covariance matrix and then a correlated x^-fit is performed. For a single elimination Jackknife, the 
covariance matrix of a ratio of correlation functions is 

A^-1 ^ 

= [^^(*^) - ^(*^)] [^^(^^•) - ^(^^•)] (21) 

k=i 



^For a pedagogical presentation of fitting see the TASI lectures by D. Toussaint [64] . 
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where Rk is the desired ratio computed with the kth sample omitted from the fuU ensemble and R is 
the ratio computed on the full ensemble. 

Fitting correlation functions to the sum of p exponential functions to extract the ground state energy 
requires fitting ranges that start at time separations from the source that are large enough so that the 
p+ 1*^ and higher excited states have negligible contributions. The determination of the minimum time 
separation that can be included in the fit is sometimes subjective. Hence a systematic uncertainty due 
to the choice of the minimum time separation in the fit is included. This uncertainty is determined by 
observing the variation of the extracted results as a function of the choice of fitting interval. The final 
uncertainties include both systematic and statistical uncertainties combined in quadrature. 



2.4 Developments in Fitting Methodology 

Fitting correlation functions to a sum of exponentials is a notoriously hard problem. However, the 
fitting is greatly simplified and more robust if only the lowest energy eigenvalue is to be extracted. 
With this restriction, only a limited amount of information about the spectral properties of the theory 
can be extracted from the computed correlation functions and, in addition, control of the systematic 
uncertainties of the extracted ground state energy eigenvalue is somewhat limited. Because the sta- 
tistical uncertainties in baryon correlation functions grow exponentially with Euclidean time at large 
times, extracting the lowest energy eigenvalue at large time separations (in order to reduce the sys- 
tematic uncertainty due to the contamination from excited states) typically results in large statistical 
uncertainties. In other words one can trade statistical uncertainty growth for systematic uncertainty 
reduction. There are two approaches in resolving this problem. One is to develop better sampling 
methods to reduce statistical uncertainties in the correlation functions. Such techniques have been 
developed in simple lattice field theories such as spin systems (cluster algorithms). However, of the 
general correlation functions computed in Lattice QCD only one approach seems viable at this point. 
This is to extract as much information from the correlation functions at short time separations where 
the statistical noise does not overwhelm the signal. In this region, multiple exponentials contribute 
to the correlation functions, and although the general multi-exponential fit problem is difficult and 
not well behaved, correlation functions can be designed so that such fits are facilitated. Variational 
analysis on symmetric positive definite matrices of correlation functions has been successfully used in 
the Lattice QCD community to extract the energy eigenvalues contributing to the correlation functions 
(the Liischer- Wolff method). These methods were originally introduced in Refs. [Mll65j, and have been 
subsequently developed |^57l |66l [671 EBj • 



2.4.1 The Prony Method 

A simple and widely used method of estimating the mass of the ground state energy contributing to 



a correlation function is the EM, as given in Eq. (16). It is conventional to define the EM from the 
logarithm of the ratio of the correlation function on adjacent time-slices. It is also possible]^ to form 
a more general EM from time-slices separated by t j > 1. For exponentially decreasing signals with 
time-independent noise, this will naturally reduce the statistical uncertainty in the EM and improve 
the extraction of energy eigenvalues as it increases the "lever-arm" of the exponential. In such a case. 



the uncertainty in Meff.(t;tj) in Eq. (16) will decrease as 1/tj. Simple correlation functions involving 
pions have time-independent uncertainties, but this is not the case for baryonic correlation functions, 
whose relative uncertainties grow exponentially with time at large times. Improvements to baryon EMs, 
and ultimately the extraction of baryon masses and the energy eigenvalues in the volume, that result 
from tj > 1 have been explored [74]. In fitting an energy to an EM (and other generalizations), either 
the Bootstrap or Jackknife procedures are used to generate the covariance matrix associated with the 



^This was suggested by K. Juge in a talk at Lattice 2008, see Ref. [69^, but may have been used earlier. 
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time-slices in the range of the fit. This covariance matrix is then used to form the x^/dof and extract 



the best estimate of the mass as well as it's uncertainty by using methods described in Sec. |2.3[ 

Based on the findings of NMR spectroscopists, the EM method has been generalized to two or more 
exponential functions [70] and is found to compare favorably [73j to the variational approach. A 



detailed study of this approach in the analysis of correlation functions with small statistical uncertainties 
was performed in Ref. [74]. The method was found to be quite stable and to provide a simple and reliable 
method of determining the lowest energy eigenvalue. This is in contrast to an analysis of low-statistics 
correlation functions for which the method is quite unstable. 



2.4.2 The Matrix-Prony Method 

In Lattice QCD, it is simple and "cost" eflFective to construct a set of correlation functions with the same 
source (creation operator) interpolating field and several interpolating fields at the sink (annihilation 
operator). Multiple sinks and one source have been used for achieving the high statistics required 
to study multi-hadron systems. In this case, it is straightforward to generalize Prony's method to 
include all correlation functions with the same source, and we call this the matrix-Prony method. This 
method leads to a further reduction in the uncertainty of the extraction of the energy eigenvalues. A 
similar approach has been briefiy discussed in Ref. [75j. While the matrix-Prony method provides a 
computational simplification in the analysis of multiple correlation functions, and provides a tool with 
which to extract multiple energy eigenvalues and eigenstates from a single source, it should not be 
considered to be more reliable than the Liischer-WolflF method [63l [65] . 

Assume there are correlation functions from which the lowest lying energy eigenvalues are 
desired. If these correlation functions are a sum of exponentials they satisfy the following recursion 
relation, 

My{T + tj)-Vy{T) = Q , (22) 
where M and V are N x N matrices and y{t) is a column vector of N components corresponding to 



the N correlation functions. Eq. (22) implies that the correlation functions are 



N 



(23) 



n=l 



where Qn and = exp{mntj) are the eigenvectors and eigenvalues of the following generalized eigenvalue 
problem 

Mq = XVq . (24) 
Given the N sets of correlation functions, the energy eigenvalues can be found by determining the 



matrices M and V that are needed in order for the signal to satisfy Eq. (22). Solving Eq. (24) leads to 



the eigenvalues = exp(m^tj) and the eigenvectors Qn that are needed to reconstruct the amplitudes 
with which each exponential enters the correlation functions. A simple solution can be constructed as 
follows. First note that 



M 



t^tw 

E 

r=t 



t+tw 



y{T + tj)y{Tf -VJ2 y{r)y{.Tf = 



T=t 



Clearly, a solution for M and V is 



M 



t-\-tw 

E 

T=t 



-^ -1 



y{r + tj)y{T[ 



V 



t+tw 



T 



T=t 



(25) 



(26) 



The method is more generally referred to as Prony's method [71 after Gaspard Riche de Prony who first constructed 
it in 1795 [72]. These techniques and other related methods are known as hnear prediction theory in the signal analysis 
community. 

^^We have applied the method for A/" = 2 and 3. 
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where these inverses exist provided that the range, tw^ is large enough to make the matrices in the 
brackets full rank {tw > — 1). Once the eigenvalues, and eigenvectors Qn are determined, the 
amplitudes, can be reconstructed using t as a normalization point. The shift parameter tj can 
be used to improve stability and understand systematic effects. The above solution is equivalent to 
determining M and V by requiring that 

t+tw 

= J] [M|/(r + tj)-y|/(r)f [M?/(r + tj)-y|/(r)] (27) 

T=t 

is minimized. 

To go beyond extracting states, a second order recursion relation can be constructed and solved. 



The minimization condition of Eq. (27), augmented to contain the second order terms in the recursion, 
can be used to determine the unknown matrices. The resulting eigenvalue problem is a second order 
nonlinear generalized eigenvalue problem which is straightforward to solve. However, tests show that 
this approach is unstable and hence such analyses are restricted to extract, at most, energies from 
correlation functions. 




Figure 2: The generalized EMP for the mass of the S using a Matrix-Prony analysis [74j. 
The correlation function was calculated on the 20^ x 128 anisotropic clover gauge field con- 
figurations with ttItt ^ 390 MeV. The inner (darker) region corresponds to the statistical 
uncertainty, while the outer (lighter) region corresponds to the statistical and fitting sys- 
tematic uncertainties combined in quadrature. The inset shows the energies of both states 
extracted with the matrix-Prony method. 



To demonstrate how this method works, results from the S baryon correlation function calculated on 
the 20^ X 128 anisotropic clover gauge field configurations with ttItj; ~ 390 MeV are presented. Figure [2] 
shows the generalized effective mass plot (EMP) for the S mass as a function of time determined 
with a A^ = 2 matrix-Prony extraction, using both the smeared-smeared (SS) and smeared-point (SP) 
correlation functions. The inset shows the second extracted state in addition to the ground state. The 
extracted value of the S mass, determined by fitting in the time interval t = 11 to t = 50, is 

Ms = 0.24097 ± 0.00025 ± 0.00003 , xVdof = 0.81 . (28) 

The energy of the dominant state in Figure [2] plateaus around time-slice t = 10, and is well-defined over 
a large interval. 
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Figure 3: The generalized EMP for the mass of the S using a Matrix-Prony analysis for a 
variety of values for tj and tw 



The matrix-Prony method is currently our preferred approach in analyzing correlation functions. 
Detailed tests of this method were presented in Ref . ^74j , which yield ground state energies that are in 
agreement with those from the other methods. The procedure for fitting parameters and determining 
their statistical uncertainty has been described in Section [2[3| Systematic uncertainties can be calculated 
by performing fits over rolling windows of time-slices within the quoted overall range and considering 
the standard deviation of the central values of those fits. This is combined in quadrature with a further 
systematic uncertainty that is generated by sampling a large range of possible values of tj and tw 
and taking the standard deviation of the central values of the resulting fits. The generalized EMP for 
the S extracted with the matrix-Prony method for a variety of values of tw and tj can be seen in 
Figure [3| One further aspect of this method, that is perhaps the most appealing, is that the correlation 
functions corresponding to the eigenstates of the matrix-Prony matrix with lowest energy eigenvalue are 
dominated by the ground-state in the lattice volume, and hence have the longest plateau in the EMP. 
The EM corresponding to the eigenstates of the octet-baryons calculated on the 20^ x 128 anisotropic 
clover gauge field configurations with m^^ ^ 390 MeV are shown in Figure [i] (compared with Figure [2] 
and Figure [3|). 



2.5 Hadronic Interactions^ the Maiani- Testa Theorem and Liischer's Method 

Extracting hadronic interactions from Lattice QCD calculations is far more complicated than the de- 
termination of the spectrum of stable particles. This is encapsulated in the Maiani- Testa theorem [76j, 
which states that S-matrix elements cannot be extracted from infinite-volume Euclidean-space Green 
functions except at kinematic thresholds. This is clearly problematic from the nuclear physics per- 
spective, as a main motivation for pursuing Lattice QCD is to be able to compute nuclear reactions 
involving multiple nucleons. Of course, it is clear from the statement of this theorem how it can be 
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Figure 4: The EMPs for the masses of the octet-baryons from the eigenstates of a 2 x 2 
matrix-Prony analysis of the SS and SP correlation functions calculated on the 20^ x 128 
anisotropic clover gauge field configurations with m^T- ^ 390 MeV. 



evaded, Euclidean-space correlation functions are calculated at finite volume to extract S-matrix el- 
ements, the formulation of which was known for decades in the context of non-relativistic quantum 
mechanics [77] and extended to quantum field theory by Liischer fTHJ [79] . The energy of two particles 
in a finite volume depends in a calculable way upon their elastic scattering amplitude and their masses 
for energies below the inelastic threshold. As a concrete example consider tt^tt^ scattering. A tt^tt^ 
correlation function in the Ai representation of the cubic group [80] (that projects onto the continuum 
s-wave state amongst others) is 

= J] J]e^P-(^-^)(7r-(t,x)7r-(t,y)7r+(0,O)7r+(0,O)) . (29) 

In relatively large lattice volumes, the energy difference between the interacting and non-interacting 
two-meson states is a small fraction of the total energy, which is dominated by the masses of the mesons. 
This energy difference can be extracted from the ratio of correlation functions, GT^+T^+{p^t)^ where 

and where the arrow denotes the large-time behavior of 0^^+^^+. For calculations performed with p = 0, 
the energy eigenvalue, E^^ and its deviation from the sum of the rest masses of the particle, AE"^, are 
related to a momentum magnitude Pn by 

/\E^ = Er, - 2m, = 2^pI + ml - 2m, . (31) 
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To obtain k cot 5{k), where S{k) is the phase shift, the square of p„ is extracted from the energy shift 
and inserted into [771 [7H1 [7Ul IHT] 



A; cot S{k) 




S{x) 



|j|<A 

E 

j 



47rA 



X 



(32) 



where k = and which is only vahd below the inelastic threshold. The regulated three-dimensional 
sum [82] extends over all triplets of integers j such that |j| < A and the limit A ^ cx) is implicit. 
Therefore, by calculating the energy-shift, AE"^, of the two particles in the finite lattice volume, the 
scattering phase-shift is determined at AE^- In the absence of interactions between the particles the 
energy eigenstates in the finite volume occur at momenta p = 2ti]/ L. It is important to re-emphasize 



that the relation in Eq. (32) is valid relativistically [78l[79]. Perhaps most important for nuclear physics 
is that this expression is valid for large and even infinite scattering lengths [82j . The only restriction is 
that the lattice volume be much larger than the range of the interaction between the hadrons, which 
for two nucleons, is set by the mass of the pion. 

For the scattering of two nucleons, the scattering length is known to be unnaturally large at the 



physical pion mass, and therefore, the relation in Eq. (32) will have to be used to extract the scattering 
parameters. For systems that are not finely-tuned, such as the 7r+7r+ system, an expansion in the 
volume can be used. In the large volume limit (L ^ \a\) the energy of the two lowest-lying continuum 
states in the Ai representation of the cubic group [80] are [THl [79] 



AEi 



4:7ra 

47r^ 
ML2 



a fa 
1 - ^^L + ^Hl 



+ 



12tan5i 



ML2 



[ l + c'itan^o + c'atan^^o + ••• ] + 0{L-^) 



(33) 



where 8^ is the s-wave phase-shift evaluated at p = 27r/L. The coefficients, which result from sums 
over the allowed momenta [79] in the finite cubic volume, are Ci = —2.837297, C2 = +6.375183, 
c; = -0.061367, and ^ = -0.354156 In addition, for a > with an attractive interaction a bound 
state exists with energy (in the large volume limit) 



11 
M 



1 + 



12 



7L 1 - 27(pcot(5)' 



where (pcot5)' = ^ pcotS evaluated at = —7^. The quantity 7 is the solution of 



7 + pcoi5\p2._ 







(34) 



(35) 



which yields the bound-state binding energy in the infinite- volume limit. As expected, the finite volume 
corrections are exponentially suppressed by the binding momentum This is consistent with the 



corrections to a single particle state where the lightest hadronic excitation is the zero-momentum two- 
particle continuum state, as opposed to a state containing an additional pion for, say, the finite volume 
corrections to the single nucleon mass. 

In the limit where L <C | (p cot (^)~'^ | , which is a useful limit to consider when systems have unnaturally- 



large scattering lengths, the solution of Eq. (32) gives the energy of the lowest-lying states: 
47r2 



A^n 



ML2 



[^1 + ^2 Lp cot 5o + • • • 



47r" 
ML2 



[ d[ + 4 Lp^oiS^ + ... ] , (36) 



^^We use the nuclear physics sign convention for the scattering length. 

^^The finite volume dependence of bound states has been explored numerically in Ref. 



16 



p cot 6 (MeV) 





/ 60 
















■yio 






-0.2 J 


-20 


0.1 


02^ 


^""^^ -60 







p cot 8 (MeV) 




p cot (MeV) 




( P/«^;r )' 



Figure 5: The functions S(x) (curved orange lines) and ]9cot(^ (straight red and blue lines) 
versus (p/m^^)^ at the physical pion mass. The lower (blue) straight line corresponds to using 
the experimentally determined ^S'l-channel nucleon-nucleon scattering length and effective 
range, while the upper (red) straight line corresponds to using a scattering length of opposite 
sign but equal magnitude. The left, center and right panels corresponds to lattices with 
spatial extent of L = 8.5 fm, 24.5 fm and 36.8 fm, respectively. The intercepts of the curves 
corresponds to the energy eigenvalues in the finite volume. 



where the coefficients are rfi = -0.095901, rfs = +0.0253716, d!^ = +0.472895, 4 = +0.0790234 and 
where pcotSo in the expression for AEq is evaluated at an energy AE = di, while pcoiS^) in the 

expression for AEi is evaluated at an energy AE = d[. The values of the are determined by 
zeros of S(x), and the expressions for AEi and AEi^ excluding AE_i^ are valid for both positive and 
negative scattering lengths. 



2.6 Bound- States Versus Scattering- States 

It is important to understand what can be extracted from finite-volume calculations. One important 
question that arises is: if a negative energy-shifted state is calculated on the lattice at finite-volume, 
does it correspond to a bound state or to a scattering state? Clearly, calculations of the same correlation 
function in multiple volumes will allow for the exponential volume dependence of a bound state to be 
distinguished from the power-law volume dependence of a scattering state. However, one can make 
an educated guess about the nature of the state by the magnitude of the energy shift. Consider a 
simple system whose scattering amplitude is dominated at low-energies by the scattering length and 
effective range, as is the case for the scattering of two nucleons. The location of the states in the lattice 



volume are determined by the solution of Eq. (32). In Figure [5| we show the graphical solution to 
Eq. (32) for two systems, one with a = +5.425 fm and r = 1.749 fm (blue line) and the other with 
a = —5.425 fm and r = 1.749 fm (red line). One finds that states with < and pcotS < (which 
occur for x = (|^) < di) likely correspond to a bound-state, while states with < and pcotS > 
{x > di) likely correspond to continuum states. However, one can imagine scattering parameters in the 
momentum expansion of pcot S that modify these rules. 

The location of the bound-state and the first continuum state for two nucleons in the ^Si channel 
(neglecting D-wave interactions and mixing) as a function of lattice volume are shown in Figure [6} In 
order for the bound-state energy to be very close to that of the deuteron binding energy, the lattice 
volumes must be very large. However, as shown in Section [6} a single calculation of the lowest energy in 
the lattice volume is not the best way to determine the deuteron binding energy. The deuteron binding 
energy will best be determined from the calculation of ground-state and continuum state energies in 
one or more lattice volumes. 
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Figure 6: The energy of the bound-state and first continuum state of two nucleons in the 
^Si channel (neglecting D-wave interactions and mixing) as a function of the spatial extent 
of the lattice (fm). 



2.7 Scattering Parameters from Wavef unctions. 

One remarkable feature of nuclear physics is that one can understand and compute (to reasonable 
accuracy) the properties and interactions of nuclei working with an energy-independent two-nucleon 
potential alone. Phenomenologically, one finds that the three-nucleon interaction, and, in the case 
of electroweak matrix elements, meson exchange currents, are required to improve agreement with 
experiment. But the fact remains that the two-nucleon potential is the dominant interaction in nuclei. 
There is a desire to construct such nucleon-nucleon potentials directly from QCD, and hence from 
Lattice QCD. Nucleon-nucleon potentials may be defined from Lattice QCD calculations in the same 
way that phenomenological potentials are determined from experimental measurements of the elastic 
scattering cross-section. A large number of Lattice QCD calculations will be performed, producing 
values for the phase shift, along with an uncertainty, over a wide range of low-energies. Potentials can 
be defined that minimize the x^/dof in a global fit to the calculated phase-shifts. At present, there is 
no practical program underway to perform such an analysis due to limited computational power. 

The utility of two particle wavefunctions and potentials in LQCD calculations was first addressed 
in a substantive way by Liischer [79j. The methodology outlined by Liischer has been used by the 
CP-PACS collaboration in early work that determined scattering parameters and phase shift of the 
TT^TT^ system from quenched calculations at relatively large pion masses with Wilson fermions [841 [85] . 
Consider an interpolating operator for n^n^ states of the form 

^(x, y; t) = ^(x, t)-f^d{^, t)u{y, t)75rf(y , t) , (37) 

where the quark-field operators q{x) may be smeared about the point x to enhance the overlap of the 
single- 7r+ interpolating operator, t^(x, t)75rf(x, t), onto the single-7r+ state. A tt^tt^ correlation function 

^ Z^^{0;k) Z^^{\^\;k) e-^^°^' tP^^{\^\;k) , (38) 

can be constructed, where R represents an element of the cubic group. The magnitude of momentum, 
k, i. defined through = 2 ^FT^, a. in Eq. El. The =ummatta over R and y projects onto the 



^^In this context, the word potential means an energy-independent potential. 
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Ai representation of the cubic-group with zero total momentum. We have performed a separation into 
the product of a sink-dependent overlap factor, Z^^{\:s.\] k)^ and a hadronic wavefunction, k). 
The spatially-dependent overlap factor Z^^{\:^\]k) is introduced in order to stress the fact that the 
overlap (onto a given eigenstate in the lattice volume) of the composite sink-operator ^(R[x] +y, y; t) is 
not the square of the overlap of the interpolating operator for the single when |x — R where R 
is a typical strong interaction length scale set by the mass of the pion (i.e. it is set by the range of the 



interaction between two tt+'s.). The wavefunction given in Eq. (38), (/)(|x|;/c) ^ Zt^t^{\:s.\] k) iIjt^t^{\:s.\] k) 
satisfies the Schrodinger equation \J9\ 

{V^ + k^)(t>{\^\;k) = jd''yUk{^,y)(t>{\y\;k) , (39) 

where Uk{'^^ y) is the Fourier transform of the modified Bethe-Salpeter kernel for the tt^tt^ interaction, 
which is generally nonlocal and energy-dependent [791 IHH [85] . An investigation of this construction 
under the assumption that [//^(x, y) is independent of energy can be found in Ref. [86j. In the limit 
where the separation between the two tt+'s is large compared to the range of the interaction (the in- 
verse pion mass), |x| i?, the normalization factor tends to constant value that is the product of the 
single- 7r+ overlap factor, Z^^{\:s.\]k) Z^^k^^ and the non-local, energy-dependent kernel tends to 
zero, Uki'x^y) 0, leaving the hadronic wavefunction, k)^ to satisfy the Helmholtz equation, 

( + /c^ ) ^/^TT-TT-dxl; fc) = 0. As the Ai representation of the cubic-group projects onto angular mo- 
mentum / = 0, 4, 6, 8, .., at low-momentum only s-wave scattering is significant and outside the range of 
the interaction the wavefunction becomes ^/^T^T^dxl; k) = Aq ^0(^1x1) + 5o n/(/c|x|). The scattering phase 
shift is simply determined from tan5o(fc) = —Bq/Aq. The CP-PACS collaboration has claimed that 
the uncertainty in the scattering length extracted from LQCD calculations of the asymptotic behavior 
of the TT^TT^ wavefunction is approximately two-thirds that of the extraction from the energy eigenvalue 
via the usual Liischer-method |8ll [85] . 

One can define an energy- and sink-dependent potential, defined at one energy (the energy of the 
two TT^'s determined in the finite volume) through an uncontrolled modification to Eq. (|39|), 



V^ + .^)^(N;AO_ ^^^1^1, _ (40) 



0(|x|;/c) 

but it is clear that, except in the case of infinitely massive hadrons [87], V/e(|x|) contains no more 
information than the phase shift at the energy determined with Liischer's method One such po- 
tential for nucleon-nucleon scattering was calculated in quenched QCD in Ref [88j but, for the reasons 
detailed above, it cannot be used as an input in nuclear calculations and it cannot be meaningfully 
compared to phenomenological nucleon-nucleon potentials An investigation of the sink dependence 

^^Alternatively, we could have simply indicated that the wavefunction depends upon the sink by k) 
4t^(|x|;A:). 

The phase-shift is recovered by solving 

( -V^+m^ y,(|x|) )(/)(|x|;fc) = e^{\^\;k) , (41) 

at the calculated value of k'^. 

The information content of the potentials derived from the analysis of hadron-hadron wavefunctions calculated with 
Lattice QCD is displayed by the following modifications to the method of Ref [88 . One could use wall-sinks for the two 

hadrons, J7^(0,t), by projecting each quark to have zero momentum. The interpolating operator has the same quantum 
numbers as that used in Ref |88j and the associated correlation function 

G(x,t) = {0\j\0,t)J{0,0)\0) , (42) 

is uniform over the spatial volume at each time-slice, and gives rise to a potential that is also uniform on each time-slice: 
Vfc(|x|) = k'^/m. This provides a clear demonstration of the fact that potentials extracted by such methods depend upon 
the choice of sink operator and are not unique predictions of QCD. 
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Figure 7: Present constraints on threshold s-wave tttt scattering. Noteworthy in the left 
panel [i93j are the red hatched ellipse from the Roy equation analysis and the grey band from 
the direct Lattice QCD calculation of the n^n^ scattering length, as discussed in the text. 
The right panel shows the n^n^ scattering length results only (the ETM result is at the top, 
the NPLQCD result is in the middle and the CGL result is at the bottom). 



of these potentials has been performed and presented in a set of lectures [89j . In this study, block spin 
transformations were applied to the quark fields used in the calculations and 0{1) modifications to the 
potentials were found out to a distance of ^ 0.5 fm. Analogous energy- and sink-dependent potentials 
have been constructed for hyperon-nucleon scattering as well [90] . 



3 Two-Body Physics 

3.1 Meson-Meson Interactions 

The low-energy scattering of pions and kaons, the pseudo-Goldstone bosons of spontaneous chiral sym- 
metry breaking, provides a perfect testing ground for Lattice QCD calculations of scattering parameters. 
There is little or no signal-to-noise problem in such calculations and therefore highly accurate Lattice 
QCD calculations can be performed with moderate resources. Moreover, the EFTs which describe the 
low-energy interactions of pions and kaons, including lattice-spacing and finite- volume effects, have been 
developed to non-trivial orders in the chiral expansion. 

The 1 = 2 pion-pion (tt^tt^) scattering length serves as a benchmark calculation with an accuracy 
that can only be aspired to at present for other systems. Furthermore, due to the chiral symmetry of 
QCD, TTTT scattering at low energies is the simplest and best-understood of the hadron-hadron scattering 
processes. The scattering lengths for tttt scattering in the s-wave are uniquely predicted at LO in x- 
PT 

m^4^^ = 0.1588 ; m^a^^^ = -0.04537 , (43) 

when is set equal to the charged pion mass. While experiments do not directly provide stringent 
constraints on the scattering lengths, a determination of s-wave tttt scattering lengths using the Roy 
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equations has reached a remarkable level of precision [92^^ ^93] : 

m^4^^ = 0.220 ±0.005 ; m^a^^^ = -0.0444 ± 0.0010 . (44) 

The Roy equations [94j use dispersion theory to relate scattering data at high energies to the scattering 
amplitude near threshold. At present, Lattice QCD can compute tttt scattering only in the 1 = 2 channel 
with precision as the / = channel contains disconnected diagrams which require large computational 
resources. It is of great interest to compare the precise Roy equation predictions with Lattice QCD 
calculations. Figure [7| summarizes theoretical and experimental constraints on the s-wave tttt scattering 
lengths [93]. It is clearly a strong-interaction process where theory has somewhat out-paced the very 
challenging experimental measurements. 

The only existing n/ = 2 + 1 Lattice QCD prediction of the / = 2 tttt scattering length involves a 
mixed-action Lattice QCD scheme of domain-wall valence quarks on a rooted staggered sea. Details of 
the lattice calculation can be found in Ref. [95]. The scattering length was computed at pion masses, 

- 290 MeV, 350 MeV, 490 MeV and 590 MeV, and at a single lattice spacing, b - 0.125 fm and 
lattice size L ^ 2.5 fm [95j. The physical value of the scattering length was obtained using two-flavor 
mixed-action x-PT which includes the eflFect of flnite lattice-spacing artifacts to 0{rri^h^) and (9(6^) [43j . 
The flnal result is: 



m^^f = -0.04330 ± 0.00042 , (45) 

where the statistical and systematic uncertainties have been combined in quadrature. Notice that 
this is a 1% calculation, but it does rely on the assumption that the rooting procedure used in the 
generation of the staggered gauge fleld conflgurations is valid, and also relies on mixed-action chiral 
perturbation theory at NLO to describe the lattice spacing artifacts. The agreement between this 
result and the Roy equation determination is a striking conflrmation of the lattice methodology, and 
a powerful demonstration of the constraining power of chiral symmetry in the meson sector. However, 
lattice calculations at one or more smaller lattice spacings are required to verify and further reflne this 
calculation. It is also desirable that the scattering amplitude be calculated with mixed-action chiral 
perturbation theory to one higher order in the lattice spacing in order to explore the convergence of 
this expansion. 

It is of great importance to have other Lattice QCD determinations of the s-wave meson-meson 
scattering lengths which use different lattice discretizations. Very recently, the ETM collaboration 
has performed a n/ = 2 calculation of the / = 2 tttt scattering length at pion masses ranging from 

^ 270 MeV to 485 MeV, at two lattice spacings (6 ^ 0.086 fm and b ^ 0.067 fm) and at two lattice 
sizes (L ^ 2.1 fm and L ^ 2.7 fm) [96]. The result extrapolated to the physical pion mass is: 



-0.04385 ± 0.00028 ± 0.00038 , (46) 



where the flrst uncertainty is statistical and the second is an estimate of systematic effects. The 
agreement among the Lattice QCD calculations and the Roy equation determination is striking. In 
Figure |8l the CP-PACS, NPLQCD and ETM Lattice QCD calculations of the 7r+7r+ scattering length 
are shown with the LO X"PT prediction subtracted. The shaded band corresponds to the NLO X"PT 
flt to the lattice calculations. 

It is interesting to consider the pion mass dependence of the meson-meson scattering lengths 
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as compared to the current algebra predictions. In Figure M (left panel) one sees that the / = 2 tttt 



The K^K^ and 7r+i^+ scattering lengths have also been computed by the NPLQCD collaboration. We refer the 
interested reader to Figure [9] and Figure 10, and to Ref. [97J for details. Exploratory calculations of ttK scattering in 



both isospin channels have been recently performed using rif = 2 + 1 improved Wilson quarks [98 , and two quenched 
studies have also been performed ^99l llQQj . 
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Figure 8: Lattice QCD calculations of the tt+tt^ scattering length with the LO x-PT pre- 
diction removed [96j. The shaded band corresponds to the NLO x-PT fit to the lattice 
calculations. The calculations with L = 2.1 fm and 2.7 fm were performed by the ETM 
collaboration [96] . 



scattering length is consistent with the current algebra result up to pion masses that are expected 
to be at the edge of the chiral regime in the two-fiavor sector. While in the two fiavor theory one 
expects fairly good convergence of the chiral expansion and, moreover, one expects that the effective 
expansion parameter is small in the channel with maximal isospin, the lattice calculations clearly imply 
a cancellation between chiral logs and counterterms (evaluated at a given scale). However, as one sees 
in Figure [9] (right panel), the same phenomenon occurs in K^K^ where the chiral expansion is governed 
by the strange quark mass and is therefore expected to be much more slowly converging. The n^K^ 
scattering length exhibits similar behavior when /llktt <^K+7r+ is plotted against /xxtt/V/k/tt as shown 
in Figure [TOl This remarkable cancellation between chiral logs and counterterms for the meson-meson 
scattering lengths is quite mysterious. 

One (very) naive interpretation of these results is that the contributions to these quantities from 
higher order terms in the chiral expansion are much smaller than one would naively anticipate, and this 
is not just a result of cancellation between terms. This leaves one with the task of understanding the 
origin of this suppression of the higher order terms. The totality of precision results may suggest that 
there is a dynamically induced length scale governing the size of contributions beyond tree-level that 
has yet to be understood. 



3.2 Meson-Baryon Interactions 

Pion-nucleon scattering has long been considered a paradigmatic process for the comparison of x~ 
PT and experiment. To this day, controversy surrounds determinations of the pion-nucleon coupling 
constant and the pion-nucleon sigma term. The K~n interaction is important for the description of kaon 
condensation in the interior of neutron stars and meson-baryon interactions are essential input 

in determining the final-state interactions of various decays that are interesting for standard model 
phenomenology (see Ref. [102j for an example). In determining baryon excited states on the lattice, it 
is clear that the energy levels that represent meson-baryon scattering on the finite- volume lattice must 
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Figure 9: mT^a^+T^+ vs. m^/ (left panel) and mkaK+K+ vs. rrik/ fk (right panel). The black 
(darker) circles and dark-brown triangles are the results of Lattice QCD calculations by the 
NPLQCD collaboration, the blue (darker) square (left panel) is from the CP-PACS collabo- 
ration, the purple (lighter) circles, red diamonds and magenta (lighter) squares (left panel) 
are from the ETM collaboration. The solid (red) lines are the current algebra predictions. 

nK^ (1=3/2) 





-0.05 



+ -0-1 

+ 

^ -0.15 



-0.2 



-0.25 



0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 




Figure 10: /xx7r<^K+7r+ vs. /xxtt/V/k/tt- The points are the results of Lattice QCD calcula- 
tions by the NPLQCD collaboration and the solid (red) line is the current algebra prediction. 



be resolved before progress can be made regarding the extraction of single-particle excitations. 

While pion-nucleon scattering is the best-studied meson-baryon process, both theoretically and 
experimentally, its determination on the lattice is computationally prohibitive since it involves anni- 



hilation diagrams , Combining the lowest-lying meson and baryon flavor octets, one can form flve 



meson-baryon elastic scattering processes that do not involve annihilation diagrams: 7r+S+, tt+S^, K^p^ 



^^In recent work, the s-wave ttN phase-shift was extracted from the negative-parity single nucleon correlation func- 
tion 
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K^n^ and K^E^ Three of these processes involve kaons and therefore are, in principle, amenable to 
an SU (3) HBx-PT analysis [103j for extrapolation. The remaining two processes involve pions inter- 
acting with hyperons and therefore can be analyzed in conjunction with the kaon processes in SU{3) 
HBx-PT, or independently using SU{2) HBx-PT. An analysis of meson-baryon scattering using the 
mixed-action technology and resources was recently performed (104j . The scattering lengths of the 
five meson-baryon processes without annihilation diagrams have been calculated to (9(m^^) in SU{3) 
HBx-PT [105. .106J. and involve low-energy constants, the Cs and the /I's, and the loop functions, 3^'s, 
which are given in Ref. jl04j . The system of processes is found to be over-constrained, and multiple 
fitting strategies are possible, as discussed in Ref. |104j . It is convenient to rewrite the x-PT formulas as 
polynomial expansions. For instance, in the case of tt^S^, the NLO and next-to-next-to-leading-order 
(NNLO) polynomial expressions are 



Pnlo 
Pnnlo 



+ S0 



47r/^a^+so 



1 + 



1 + 



Ms 



= 1 - Coim^ 

777/77- 



(47) 



The left-hand sides of these equations are given entirely in terms of lattice-determined quantities, while 
the right-hand side provides a convenient polynomial fitting function. Figure [TT] shows the results of 
the Lattice QCD calculation, along with the fits. The shift of the value of F from NLO to NNLO is 
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Figure 11: Fnlo and Fnnlo versus the meson masses for two of the meson-baryon processes in SU(3) 
HBx-PT. The line at F = 1 is the LO curve, and dotted line is the physical meson mass. 



dependent on the renormalization scale /x, and therefore with the choice /x = one would expect this 
shift to be perturbative if the expansion is converging. The large shifts in F from NLO to NNLO are 
indicative of large loop corrections. The low-energy constants (LECs) fit to the results of the lattice 
calculations are tabulated in Ref. |104j . While the NNLO LECs hi and hus appear to be of natural 
size, the NLO LECs Cq and Cqi are unnaturally large. The extrapolated values of the five scattering 
lengths are given in Table [T] While the 7r+S+ and tt+S^ scattering lengths appear to be perturbative, 
the extrapolated kaon-baryon scattering lengths at NNLO deviate by at least 100% from the LO values. 
The seemingly inescapable conclusion is that the kaon-baryon scattering lengths are unstable against 
chiral corrections in the three-fiavor chiral expansion, over the explored range of light-quark masses. 

Given the poor convergence found in the three-fiavor chiral expansion due to the large loop correc- 
tions, it is natural to consider the two-fiavor theory with the strange quark integrated out. In this way, 

has the same quantum numbers as tt+S^. 
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Table 1: SU{3) HBx-PT extrapolated scattering lengths. The first uncertainty is statistical, 
and the second is the statistical and systematic uncertainty added in quadrature. "NLO 
(NNLO fit)" indicates that the Ci and Cqi from the NNLO fit to 7r+S+ and tt+S^ have been 
used. 



Ouantitv 


LO (fm) 


NLO fit (fm) 


NLO (NNLO fit) (fm) 


NNLO (fm) 




-0.2294 


-0.208(01)(03) 


-0.117(06)(08) 


-0.197(06)(08) 


ans 


-0.1158 


-0.105(01)(04) 


0.004(05)(11) 


-0.096(05)(12) 




-0.3971 


-0.311(18)(44) 


0.292(35) (48) 


-0.154(51)(63) 


0-Kn 


-0.1986 


-0.143(10)(27) 


0.531(28)(68) 


0.128(42)(87) 




-0.4406 


-0.331(12)(31) 


0.324(39) (54) 


-0.127(57)(70) 



ttS and ttS may be analyzed in an expansion in m^j^. To (9(m^) in the two-fiavor chiral expansion, one 
has [TOTj 
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(48) 



47r 777.7,- + 

where ^4 is the LEG which governs the pion mass dependence of /^^ j92] . Note that the chiral logs have 
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Figure 12: a versus the pion mass for 7r+S+ and tt+S^. The diagonal line is the LO curve, and the 
dotted line is the physical pion mass. The innermost error bar is the statistical uncertainty and the 
outermost error bar is the statistical and systematic uncertainty added in quadrature. The filled bands 
correspond to the fits to the LECs in the SU(2) case at NNLO. 

canceled, and in this form, valid to order in the chiral expansion, the scattering lengths have a 
simple polynomial dependence on m^^ |107j . Exploring the full 95% confidence interval error ellipse in 
the h-C plane yields pion-hyperon scattering lengths, extrapolated to the physical pion mass, of 



-0.197 ±0.017 fm 



-0.096 ±0.017 fm 



(49) 



and the scattering length versus the pion mass is shown in Figure 12 
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The bar denotes the scattering length rescaled by a kinematical factor [T 
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The HBx-PT analyses of meson-baryon scattering support a general observation about convergence 
in the three-flavor chiral expansion. As the pion masses considered in the lattice calculation are com- 
parable to the physical kaon mass, the distinct convergence patterns of the two- and three-flavor chiral 
expansions are suggestive that the breakdown in the three-flavor case is not due to the relative large- 
ness of the strange-quark mass as compared to the light quark masses, but rather due to some other 
enhancement in the coeflicients of the loop contributions, possibly related to a scaling with powers of 
n/, the number of flavors. 

3.3 NN, YN and YY Interactions 

Perhaps the most studied and best understood of the two-hadron processes are proton-proton and 
proton-neutron scattering. In the S-wave, only two combinations of spin and isospin are possible, a 
spin-triplet isosinglet np {^Si) and a spin-singlet isotriplet pp (^Sq). At the physical pion mass, the 
scattering lengths in these channels are unnaturally large and the ^Si —^Di coupled-channel contains 
a shallow bound state, the deuteron, with a binding energy of ^ 2.22 MeV. These large scattering 
lengths and the shallow bound state are described, in EFT, by the coeflicient of the momentum- 
independent four-nucleon operator having a non-trivial flxed-point in its renormalization group flow 
for the physical light-quark masses. An interesting line of investigation is the study of the scattering 
lengths as a function of the quark masses to ascertain the sensitivity of this flne-tuning to the QCD 
parameters |108[ I109[ IllOj . The flne tuning is not expected to persist away from the physical masses 
and we expect present day (unphysical) Lattice QCD calculations to yield scattering lengths that are 
natural-sized. More generally, it is interesting to determine how the structure of nuclei depends upon 
the fundamental constants of nature. In particular, we expect that any nuclear observable is essentially 
a function of only flve constants, the length scale of the strong interactions, Aqcd^ the quark masses, 
rUu^ rrtd and m^, and the electromagnetic coupling ^ 

The flrst study of baryon-baryon scattering with Lattice QCD was performed more than a decade 
ago by Fukugita et al I112j . This calculation was quenched and at relatively large pion masses, 
m^^^ 550 MeV. Since this time, the dependence of the NN scattering lengths upon the light-quark 
masses has been determined to various non-trivial orders in the EFT expansion [108^ lOQilllOj . which is 
estimated to be valid up to ttIt^ ~ 350 MeV. Therefore physical predictions of NN scattering parameters 
becomes possible with Lattice QCD calculations that are performed with pion masses less than ^ 
350 MeV. 

The NPLQCD collaboration performed the flrst n/ = 2 + 1 QCD calculations of nucleon-nucleon in- 
teractions [113] and hyperon-nucleon [114J interactions at low-energies but with unphysical pion masses, 
and the nucleon-nucleon scattering lengths were found to be of natural size. The flne-tunings at the 
physical values of the light-quark masses indicates that Lattice QCD calculations with quark masses 
much closer to the physical values (than today) are needed to extrapolate to the experimental values. 
The results of the Lattice QCD calculation at the lightest pion mass and the experimentally-determined 
scattering lengths at the physical value of the pion mass were used to constrain the chiral dependence 
of the scattering lengths from ttIt^ ~ 350 MeV down to the chiral limit [113j. However, these results 
suggest various possible scenarios toward the chiral limit which can only be resolved by way of Lat- 
tice QCD calculations at lighter pion masses. In contrast, very little is known about the interactions 
between nucleons and hyperons from experiment, and Lattice QCD calculations can provide the best 
determinations of the corresponding scattering parameters and hence determine the role of hyperons in 
neutron stars. 

The energy eigenstates in the flnite lattice volume are classifled by their global quantum numbers, 
baryon number, isospin, third component of isospin, strangeness, total momentum, and behavior un- 

^^In the low-energy theory, the dependence on top, bottom and charm quark masses is encapsulated in Aqcd- 
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Table 2: Baryon-baryon channels calculated in Ref. |117j . 



Channel 



PP {'So) 
np {^Si) 
nA {'So) 
nA {^Si) 
nS- {'So) 

S-S- {'So) 

AA {'So) 
S-S- {'So) 



+1 






-1 






-1 
-1 
-1 
-1 

-2 
-2 
-4 



der hyper-cubic transformations. Six quark operators that are simple products of three-quark baryon 
operators are (generally) used as sources for the baryon-baryon correlation functions, but this is not a 
requirement. As a consequence, the baryon content of the interpolating operator is used to define the 
operator, e.g. nA(^S'i), but this operator will, in principle, couple to all states in the volume with the 
quantum numbers 5 = 2, / = |, 1^ = — |, s = —1, and ^^^'Lj '^Si + where the ellipses denote 
states with higher total angular momentum that also project onto the Ai irreducible representation of 



the cubic-group Correlation functions for the nine baryon-baryon channels shown in Table [2] have 
been calculated, using both SS and SP correlators, as described in Ref. If the calculations were 

performed on gauge field configurations of infinite extent in the time-direction, so that only forward 
propagation could occur, some of the channels in Table [2] could be analyzed by considering contributions 
from a single scattering channel, e.g. A^A^, nS~, for which a single, well-separated ground 

state with these quantum numbers is expected. However, other channels may require a multi-channel 
analysis, e.g. nA, AA. The nA source will produce low-lying states in the lattice volume that are 
predominately linear combinations of the nA, nTP and pS" two-baryon states. The A A source will 
produce low-lying states in the lattice volume that are predominately linear combinations of the AA, 
S^'^S^'^, and A^S two-baryon states. 

There are a number of ways to perform the statistical analysis of the correlation functions and 
determine the associated uncertainties. One way is to use the Jackknife method to determine the 



uncertainty in pcotS directly. However, this is complicated by the fact that S{x) in Eq. (32) is a 
singular function. An alternate method is to determine the value of and its associated uncertainty 



from the two-baryon energy splitting from Eq. (31), and then to propagate the central value and la 



uncertainties through Eq. (31) to determine pcotS. 

An extensive exploration of the impact of high-statistics on one-, two- and three-baryon correlation 
functions on one ensemble of anisotropic clover gauge configurations generated by the Hadron Spectrum 
Collaboration was undertaken in Refs [TU 11161 1117 ]. A total of ^ 440,000 sets of calculations were 
performed using 1200 gauge configurations of size 20^ x 128 with an anisotropy parameter ^ = bg/bt = 



The spatial dimensions of the gauge field configurations that have been used to date for such calculations are identical, 
and as such the eigenstates of the QCD Hamiltonian can be classified with respect to their transformation properties under 
cubic transformations, i^(3), a subgroup of the group of continuous three-dimensional rotations, 0(3). The two-baryon 
states that are calculated in this work all belong to the Af representation of i^(3), corresponding to states with angular 
momentum L = 0, 4, 6, . . . . 
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3.5, a spatial lattice spacing of bg = 0.1227 ± 0.0008 fm, and pion mass of ^ 390 MeV. The ground 
state baryon masses (in lattice units) were extracted with uncertainties that are at or below the ^ 0.2%- 
level. Figure 13 shows the effective |kp plot for both the proton-proton and neutron-proton channels. 
Both channels exhibit plateaus in |kp. While the plateau in the proton-proton channel differs from zero 




Figure 13: The left (right) panel is the effective |kp plot for the proton-proton {^Sq) (neutron- 
proton (^S'l)) channel and the fit to the plateau for the calculations on anisotropic clover 
gauge configurations in Ref. [117J. 



at the ^ 1-cr level, the plateaus in both channels are consistent with zero. We conclude that at this value 
of the pion mass, the interactions between nucleons produce a small scattering length in both channels 
compared to the naive estimate of ^ 0.5 fm. A summary of all Lattice QCD calculations of NN 
scattering is shown in Figure [T4j The results calculated on the anisotropic clover gauge configurations 
are consistent with those obtained with mixed- action Lattice QCD jll3j . It is interesting to note that 
the results of quenched calculations |118j yield scattering lengths that are consistent within uncertainties 
with the fully-dynamical n/ = 2 + 1 values. Figure 15 shows a summary of the high-statistics results 
for baryon-baryon interactions obtained on the anisotropic clover gauge configurations. The AA channel 
is found to be negatively shifted in energy which may signal that the lowest-state is in fact a bound- 
state, but calculations in a larger volume are required before more definitive conclusions can be drawn. 
The calculations constitute an order of magnitude jump forward in the volume of output for this type 
of calculation. The two additional lattice volumes in which we are presently performing calculations 
at this pion mass will allow for a systematic exploration of the volume dependence of the scattering 
amplitude for all of the two-hadron systems. In principle, this will enable a separation of scattering 
and bound states. However, calculations at this pion mass should be viewed to be only a proof of 
principle and a testing ground for new analysis techniques. In order for lattice calculations to provide 
meaningful constraints on interactions at the physical pion mass through chiral extrapolation, the pion 
mass must be substantially reduced down toward the physical value {m^ ^ 140 MeV) while maintaining 
the integrity of the calculation (i.e. small enough lattice spacings and large enough volumes). 
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Figure 14: A compilation of the scattering lengths for NN scattering in the ^S^ (left panel) 
and ^aSi (right panel) calculated with Lattice QCD and with quenched Lattice QCD. The 
vertical lines correspond to the physical pion mass. 




Figure 15: Baryon-baryon interactions extracted from calculations at m^^ ^ 390 MeV and 
in a lattice volume of 20^ x 128 |117j . The top-most point of the plot-legend corresponds to 
the left-most point on the plot, and the bottom- most point of the plot-legend corresponds 
to the right-most point on the plot. The other points are ordered accordingly. 



4 Few-Body Physics 

Computing resources and lattice algorithms have reached a stage where it is now possible to consider 
more complicated hadronic observables such as those in the baryon number, B > 2 sector — the domain 
of nuclei. Here there are many observables about which very little (or nothing) is known experimentally 
or theoretically, but which are phenomenologically important to nuclear structure and interactions and 
to nuclear astrophysics. Systems containing more than two mesons are of phenomenological interest in 
a number of areas from heavy ion collisions at RHIC, to the equation of state of neutron stars. In the 
last few years, there has been a concerted effort to understand and develop the techniques needed for 
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studying systems of a few or many hadrons. 



4.1 Multi-Meson Interactions 

As discussed previously, it has been known for many years how to exploit the volume dependence of 
the eigen-energies of two-hadron systems to extract infinite volume scattering phase shifts [78] provided 
that the effective range of the interaction, r is small compared to the spatial extent of the lattice volume 
(since r ~ for most interactions, this constraint is m^r L ^ 1). In recent works, this has been 
extended to systems involving n > 2 bosons |119l 11201 1121 J and n = 3 fermions |122j in the situation 
where the relevant scattering length, a, is small compared to the spatial extent of the lattice. By 
performing a perturbative EFT calculation in a finite volume, the ground state energy of a system of 
n bosons was computed in Refs. [ 1191 112()[ 1121] . The shift in energy of n bosons of mass M from the 
non-interacting system is 
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where the parameter a is related to the scattering length, a, and the effective range, r, by 
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(50) 



(51) 



The geometric constants, X, J^, /C, >C, 7o,i, that enter into Eq. (50) are defined in Ref. (121] and 
^Cm are the binomial coefficients. The three-body contribution to the energy-shift given in Eq. (50) is 
represented by the parameter T]^ (see Ref. |121j ). 

Lattice QCD calculations of these energy shifts allow for an extraction of the parameters a and 
7^3 . To determine the energy shifts, the multi-meson correlation functions (specifying to the multi-pion 
system) 



C^it) OC (^(^7r-i^,t)j |^7r+(0,0)j ^ , 



(52) 



are calculated. On a lattice of infinite temporal extent, the combination 



Gn{t) 



_ Cn{t) t- 



(53) 



^^Effects of temporal (anti-)periodicity are discussed in Ref. [12 
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allows for an extraction of the ground-state energy shift, AE"^, which can then be used as input into 



Eq. (50) to extract the scattering and interaction parameters. To compute the (n!)^ Wick contractions 



in Eq. (52), the correlation function can be written as 



^/3i/32../3 



(54) 



where aS(x, t; 0, 0) is a light-quark propagator. The object (block) 11 is a x bosonic time-dependent 
matrix where = 12 (A^ = A^^- x Nc with A^^- = 4-spin and Nc = 3-color), and rj^ is a twelve component 
Grassmann variable. Further simplifications are possible resulting in the correlation functions being 
written in terms of traces of powers of 11. As an example, the contractions for the 3-7r^ system give 



C3(t) OC trc,s [n]' - 3 trc,s [n^] trc,s [H] + 2 trc,s [n'] , 



(55) 



where the traces, tr^s, are over color and spin indices. Contractions for n < A^ tt+'s are given explicitly 
in Ref. [mi- 



4.2 Two- and Three- Body Interactions 

The NPLQCD collaboration has performed mixed- action Lattice QCD calculations of the n < 12 pion 
and kaon correlation functions |123i 11241 I125j . In order to correctly calculate these correlators for 
large n, very high numerical precision is necessary and the ARPREC arbitrary precision numerical 
library [126j was used. By performing a correlated fit to the eflFective energy diflFerences extracted 
from these calculations, the two- and three-body interactions were determined |123l 11241 1125j . and the 
two-body interactions agree with those extracted from the two-body sector alone [95j. The resulting 
three-body interactions are displayed in Figure 16, The 3-7r+ interaction is found to be repulsive with 
a magnitude consistent with the expectation from naive dimensional analysis (tree-level x-PT). In 
contrast, the three interaction is consistent with zero within somewhat larger uncertainties. 
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Figure 16: Three pion (left-panel) and kaon (right-panel) interactions determined from the 
MILC coarse (lighter, blue) and fine (darker, magenta) lattices plotted as a function of the 
dimensionless ratios m.,^/ f.,^ and ttikI fx respectively. 
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4.3 Mixed Species Meson Systems 



In the case of a mixed system comprising n mesons of one type and m mesons of a second type, the 
results summarized in Section 4.1 have been generahzed by Smigielski and Wasem |127j . working to 
0{1/L^). The energy shift of the interacting system from the free system depends on three two-body 
interaction parameters and four three-body interaction parameters. As might be expected, the fuU 
form of the energy shift is cumbersome and we refer the interested reader to the original paper. These 
expressions are currently being used to analyze LQCD calculations of the energies of systems of pions 
and kaons ^28j. Correlation functions are calculated in a way that is similar to the way that the 
multi-pion correlation functions are calculated. 



4.4 Contractions for Large (N > 12) Systems of Mesons 

In order to go beyond n = 12 meson systems, or to study systems of many different types of mesons, it 
is necessary to make use of more efficient methods of performing the contractions of quark fields. While 



better than the naive factorial construction (which scales as n!^), the construction used in Section 4.1 
scales poorly to large numbers of mesons, behaving at best as n!^/^ (provided the matrix 11 is generalized 
to give a nonzero result). In Ref. |129j . a recursive method for performing these contractions was 
developed that allows the extension of the study of meson systems to larger n and also greatly simplifies 
the contractions required for systems of many different species of mesons. We will summarize the 
construction by first considering the recursive approach to the contractions for a single species of meson 
before reporting the general case. 

By rescaling the correlation functions of the single-species, single-source system considered previously 

as 

Cn.^{t) = (-r n! {R^) , (56) 

(the angle brackets denote a trace over spin and color indices) it is straightforward to show that the 
objects, Rn satisfy an ascending recursion 

Rn+i = {Rn)Il - nRnll. (57) 



with the initial condition that i?i = 11 as defined in Eq. (54). To see how this works we explicitly 
construct the first few terms: 

- i?i n = ( n ) n - 

- 2 R2II = {ufu- (n2)n-2(n)n2 + 2n^ 
3 ( ) ( n ) - 2 ( ) , (58) 



i?2 = 


{Ri 


)n 


R2) = 


(n; 




R3 = 


{R2 


)n 


R3) = 


(n; 


)3 - 



in agreement with Eq. (55). Descending recursions also exist. 



A correlation function for a system composed of riij mesons of the i^^ species from the j^^ source at 
(yj5 0), where < i < k and < j < m, is of the form 

Cn{t) = ( f E M^^^) j - f E 

nil / \ nim / \ nki / \ rifcn 

4(yi,0)) ... 4(ym,0) ... 4(yi,0) ... 4(y^,0) ), (59) 



32 



where Mi = riij is the total number of mesons of species z, and the subscript in Cn{t) labels the 
number of each species from each source, 



n 



/ nil ni2 ... nim \ 
\ riki nk2 ... rikm J 



(60) 



The w4i(y,t) denotes a quark-level interpolating operator ^^(x, t) = q^{-x^t) 75 u{^^t)^ and it is 
straightforward to show that 



c„(t) 



H-^ ^(r„) ,(61) 



where the vj are m x A^-component Grassmann variables, and the Pij are N x N dimensional matrices, 
where N = m x which are generalizations of the 11 defined in Eq. (54) with an additional species 
index, i. They are defined as 
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(62) 



The ( ^4^ are N x N dimensional matrices, one for each fiavor, and pair of source indices, a and 
b. Af = J2i -^i is the total number of mesons in the system, with Af < N. The defined implicitly in 
Eq. (jeTl) satisfy the recursion relation 



where 



EE 

i=i j=i 



T, 



n+lr 



iij — 



/ 

V 



(63) 



0\ 



... : 

••• oy 



(64) 



and where the non-zero value is in the (i, j)*^ entry Defining Uj = riij to be the number of mesons 
from the j^^ source, it is clear that the correlation function vanishes when Uj > N for any source j. 

These recursion relations (and those for the simpler systems also studied in Ref. ^129j ) allow for the 
calculation of arbitrarily large systems of mesons. As formulated above, the systems are restricted to 



contain quarks of one fiavor but anti-quarks of any number of fiavors ^ or vice versa. Importantly the 
above algorithm scales linearly with the total number of mesons in the system. 



25 



This restriction is not necessary, but relaxing it results in much more complex sets of recursions. 
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4.5 Pion and Kaon Condensation 



The ground state of a generic system of many bosons with repulsive interactions is a Bose-condensed 
phase. The QCD systems of pions or kaons discussed above form a Bose-Einstein condensate of fixed 
third-component of isospin, and strangeness, s. It is of significant theoretical and phenomenolog- 
ical interest to investigate the properties of such systems. Theoretical efforts have used LO x-PT to 
investigate the phase diagram at low chemical potential [130j and it is important to assess the extent to 
which these results agree with QCD. In neutron stars, it is possible that it is energetically favorable to 
(partially-) neutralize the system with a condensate of K~ mesons, and an isospin chemical potential 
resulting from the isospin asymmetry of the colliding nuclei may play an important role in the study 
of the QCD phase diagram at RHIC. Numerical calculations provide a probe of the dependence of the 
energy on the pion or kaon density, and thereby allow for an extraction of the chemical potential via a 



finite difference. The results from mixed-action calculations are shown in Figure 17 and Figure 18, Also 



shown are the predictions of tree-level X"PT, which are in remarkably good agreement. This is encour- 
aging for studies of kaon condensation in neutron stars where, typically, tree level X"PT interactions 
amongst kaons, and between kaons and baryons [101], are assumed. 




Figure 17: The dependence of the isospin chemical potential on the pion density. The curves 
correspond to the predictions of tree level x-PT (dashed) |13Qj . the energy shift of Eq. (50) 
(solid) and without the three-body interaction (dotted). 



4.6 Hadronic Medium Effects 

The hadronic medium formed in the multi-meson calculations described above is naturally expected to 
infiuence the properties of other hadrons that interact with it. A phenomenologically relevant example 
is the passage of a J/^^ meson through medium as J suppression is taken to be a key signal for the 
formation of a new state of matter at high temperatures, such as the quark-gluon plasma jl31j . In 
Ref. jl32j . a simplified version of this scenario was considered, and the screening of the static quark 
potential (the potential between an infinitely massive quark-anti-quark pair) by a Bose-Einstein con- 
densate of TT+'s was computed. In order to extract the shift in the static potential caused by the presence 
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Figure 18: The dependence of the strangeness chemical potential on the kaon density. The 
curves correspond to the predictions of tree level chiral perturbation theory (dashed) [130j, 



the energy shift of Eq. (50) (solid) and without the three-body interaction (dotted). 



of the medium, the following ratio of correlation functions was calculated 



(65) 



where Cn is defined in Eq. (52) and 

CnMR^t-^tw.t) = ^0|[ J^7r-(x,t)7r+(0,tj]" W(y + r,t;y,t, 

X y,\r\=R 



(66) 



and W (y, to; y + t) is the Wilson-loop operator formed from products of gauge links joining the 
vertices at (y^to)^ (y + r,to), (y + '^it) and (y,t). At large Euclidean time, this ratio falls exponentially 
with a scale 6V{R^n) that is the shift in static potential for separation R = |r| and isospin density 
p = PqU (where po = ^ 0.064 fm~^). Further details are discussed in Ref. |132j . 

As expected from the weak nature of pion interactions with isoscalar objects, the overall screening 
effect is small. In the region where the force between the static quark and anti-quark is constant, 
F ^ 1 GeV fm~^, the shift is 0{MeV fm~^). The shift in the force was found to vary linearly with 
the isospin density of the hadronic system as shown in Figure This is consistent with a simple 
interpretation in terms of a dielectric in the volume of the color flux tube between the quark-anti-quark 
pair. Further calculations extending this study to dynamical charm quarks are underway. To directly 
connect with phenomenology, the more diflicult problem of flnite baryon density must be confronted. 



4.7 Three- Baryon Systems 

The flrst signiflcant steps towards the calculation of the properties of nuclei directly from QCD were 
taken by the NPLQCD coUaboration and the PACS-CS coUaboration [I33] during 2009. The 
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Figure 19: The in-medium contribution to the radial QQ force, 5F[R^ n), at a pion number 
density of po and 5po-The inner uncertainty associated with each point is statistical, while 
the outer is the statistical and systematic uncertainties combined in quadrature. 



NPLQCD collaboration performed an/ = 2 + 1 calculation of the correlation function with the quantum 
numbers of the strangeness s = —4, baryon number 5 = 3 system, which is labeled as for 
convenience, and also of the correlation function with the quantum numbers of the triton (or ^He), 
at a pion mass of ttIt^ ^ 390 MeV [116]. Further, the PACS-CS collaboration performed a quenched 
calculation of the system with the quantum numbers of the triton and of the a-particle at a pion mass 
of - 800 MeV [T33l . 



The NPLQCD collaboration calculation used the product of single baryon sources and sinks for the 
source and sink of the three-baryon correlation functions. In the nS^S^-channel they were the product 
of two S and one nucleon source or sink, with the spin quantum numbers constructed to produce a state 
with total J = and third component of isospin 1^ = \. In the triton channel, the product of three 
nucleon interpolating operators were used to construct an operator with J \ and with total isospin 



1 

2' 



The correlation function in the nS^S^-channel, shown in Figure 20, yields an energy-splitting 



that is consistent with zero within the uncertainties of the calculation, as shown in Figure [21 



8E^ 



4.6 ±5.0 ±7.9 ±4.2 MeV , xV^of 



2.0 



(67) 



It is very encouraging that the uncertainty in the energy-shift per baryon is ^ 3 MeV, which is smaller 
than the binding-energy per nucleon in typical nuclei, B ^ 8 MeV, and not significantly larger than 
the binding-energy per nucleon in the deuteron or triton at the physical values of the light-quark 
masses. The single energy-level fit to the EMP in Figure 21 has a x^/dof = 2.0, indicating that 
there may be additional structure in the correlation function. Including a second energy-level shifted 
by AE ^ —0.004 lattice units might provide a better description of the EMP, and would be consistent 
with the lower-energy state, S^AA, that is expected to contribute to the four low-lying eigenstates in 
the lattice- volume. However, enhanced statistics are required to determine if this is, in fact, the case. 

At present, unlike the situation in multi-meson systems [1251 11241 1123j . the analytical tools are not 
in place to use the above energy shift and those of the associated two-baryon systems to extract the 
parameters describing the relevant two- and three-body interactions. While the volume dependence of 
the simplest three-fermion systems has been studied in Ref. jl22j . the mixing we expect between four 
closely spaced states significantly complicates the situation. 

The NPLQCD collaboration focused on the state(s) that couples to the E^E^n interpolating-operator 
simply due to limited computational resources and the expectation that the large strange quark content 
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Figure 20: The EMPs for the nucleon, S, nn, E^E^ and nS^S^ systems obtained by the 
NPLQCD coUaboration [116J. 



would lead to a clean signal. In the past it has been the case that gauge field generation has required the 
majority of LQCD resources, but this is no longer true for precise calculations of baryonic observables. 
The resources that are presently required to perform the large number of calculations required for 
nuclear systems is significantly greater than that required for gauge field generation. This situation 
will improve as more eflFort is put into algorithmic improvements for contractions, in the same way that 
the use of defiation and other techniques have greatly reduced the resources required for propagator 
generation. Work in this direction is in progress. Given the observed behavior of the signal-to-noise 
ratio, which will be discussed in Section [5| identification of the ground state in systems of four and five 
baryons is anticipated in the near future. 

The number of contractions that must be performed in order to calculate the correlation function in 
the triton channel is significantly greater than the number in the S^S^n, and therefore a smaller number 
of contractions was possible with the computational resources available to NPLQCD. The results of the 
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Figure 21: The EMP of the difference between the energy of the tiEPEP channel and its 
constituent hadrons [116]. The right panel is an enlargement of part of the left panel. 
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Figure 22: The EMP for the correlation function with the quantum numbers of the triton 
(left panel), and for the difference between the triton and its constituents (right panel) [116j. 



calculations performed by the NPLQCD collaboration of the correlation function in the channel with 



the quantum numbers of the triton are shown in Figure [22j A clear plateau is found for the ground-state 
of the system in the lattice volume in the window of time-slices where the signal-to-noise ratio is not 
degrading exponentially. However, the uncertainty remains too large to determine if the ground state 
corresponds to a positively or negative shifted state, 

5Epr^n = 40 ± 21 ± 38 MeV , xV^of = 0.97 , (68) 



corresponding to ^ 13 MeV-per-baryon. However, the fact that a plateau in the triton channel has 
been observed at a relatively light pion mass {171^^ ~ 390 MeV) is a very encouraging step forward. 

Unfortunately, none of the calculations that have been performed to date for systems of three or 
more baryons have managed to calculate energy-shifts that exceed 5-a deviations from zero, and, as 
such, compelling calculations are yet to be performed. However, such calculations are expected in the 
near future. 

At the physical pion mass, the expected energy eigenvalues of two nucleons in a lattice volume with 
tUt^L ^ 1 have been determined j82]. Recently, there have been efforts undertaken to determine the 
expected energy eigenvalues of systems composed of three nucleons [ 1221 11341 11351 11361 1137j in cubic 
volumes with periodic boundary conditions using EFT. It is found that the triton binding energy is 
less sensitive to a finite-volume than the deuteron binding energy, which can be understood in terms 
of the spatial extent of the respective wavef unctions. These works are the beginning of a series of 
theoretical calculations that should be performed to provide a guide in determining the Lattice QCD 
calculations that should be performed to best extract the parameters in the EFT which will then allow 
for a description of more complex processes of relevance to nuclear physics. 



4.8 Four-Baryon Interaction 

In late 2009, the PACS-CS collaboration performed the first quenched calculation of a four-baryon 
correlation function |133j in the channel that will contain the a-particle (^He nucleus) when performed 
at the physical pion mass. The pion mass used in these calculations is large, m^^ ^ 800 MeV, and sea 
quark effects are ignored. Nonetheless this is a very important step towards calculating the properties 



and interactions of nuclei. The PACS-CS collaboration results are shown in Figure 23 One hopes to 
see significantly improved statistics in the near future for such an important quantity. 
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Figure 23: The quenched results for the binding energies (in lattice units) obtained by the 
PACS-CS collaboration in the triton channel and the channel with the quantum numbers of 
the a-particle [133j. The pion mass in these calculations is 171^^ ^ 800 MeV. 

5 The Signal-to-Noise Ratio in Baryon Correlation Functions 

Until recently it had been lore that the signal-to-noise ratio in baryon correlation functions degrades ex- 
ponentially with the baryon number. The correlation functions exhibit this behavior at large times |138j . 
and if this were the only part of the correlation functions that could be calculated with precision, the 
calculation of quantities of importance to nuclear physics would require exponentially more computing 
resources than those of importance to particle physics. However, during the last year it has been realized 
that such exponential degradation of the signal-to-noise ratio is absent at intermediate times over an in- 
terval that is dictated by the structure of the source and sink of the correlation function |115[ 11161 1117j . 
A consequence of this behavior is that it is possible to extract information about multi-baryon systems 
where the only additional computational resources required are those to perform the contractions. 

In the case of a source that has the quantum numbers of a single positive parity nucleon, the 
correlation function on an ensemble of gauge field configurations with infinite temporal extent has the 
form [138J 

{ON{t)) = 5] rf (0| A^-(x,t)]V^(O,0) |0) ^ Z^e-^-^ , (69) 

X 

where A^^(x, t) is an interpolating field (composed of three quark operators) that has non- vanishing 
overlap with the nucleon, r+ is a positive energy projector, and the angle brackets indicate statistical 
averaging over calculations on an ensemble of configurations. The variance of this correlation function 
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is given by 

= {0\ N'^{^,t)N\y,t)N\O,0)N\o,0) \0) -{9^,(1))' 



Z^^e-^^'^'-Zle-^''-' + Z^^e-^^-' + ... Z^^ e'^^-' , (70) 



where all interaction energies have been neglected, and N is the number of (independent) calculations. 
At large times, the noise-to-signal ratio has the form, as argued by Lepage |138j . 

More generally, for a system of A nucleons, the noise-to-signal ratio behaves as 

^ 1 ,A(Miv-|m^)f j<y2) 



at large times. 

The various "Z-factors", such as Z^t^^ depend upon the details of the sources and sinks interpolators 
that are used. For the calculations performed by the NPLQCD collaboration, the projection onto zero- 
momentum final state nucleons, introduces a 1 / VVolume suppression of the amplitudes of the various 
components (except for A^A^) in addition to color and spin rearrangement suppressions that exists 
independent of the spatial structure of the source. As a consequence, an interval of time slices exists at 
short times (the "Golden Window" ) in which the variance of the correlation function is dominated by 
the terms in Eq. (71) that behave as ^ g-2Mivt^ ^^iis window, the signal-to-noise ratio of the single 



baryon correlation function is independent of time. Further, the signal-to-noise ratio does not degrade 
exponentially faster in multi-baryon correlation functions than in single-baryon correlation functions in 
the "Golden Window". 

The finite temporal extent introduces backward propagating states (thermal states) into the corre- 
lation functions which lead to exponentially worse signal-to-noise ratios at large times |115i 11161 1117j . 
These contributions are suppressed by at least exp(m7T-T), however, they can cause complications. We 
note that the impact of these states can be mitigated by working at larger temporal extents and expo- 
nentially large computational resources are not required. 

With the high statistics calculations that have been performed, the behavior of the signal-to-noise 
ratio has been carefully examined, and it was found to be useful to form the effective noise-to-signal 
plot [115J. On each time slice, the quantity 

is formed, from which the energy governing the exponential behavior (the signal-to-noise energy-scale) 
can be extracted via 

E,it;t,) = ilogf^) ^ (74) 



tj s(t) 

For a correlation function that is dominated by a single state with a corresponding variance correlation 
function dominated by a single energy scale, the quantity Es{t]tj) will be independent of both t and 

The signal-to-noise ratio in the one- and two-nucleon sector has the simplest structure as only up 
and down quarks appear in the interpolating operators. In the single nucleon sector, it is expected that 
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Figure 24: The energy scale of the signal-to-noise ratio in the nucleon (left panel) and 
proton-proton (right panel) correlation functions, defined in Eq. (74), obtained on 20^ x 128 



anisotropic clover gauge field configurations [ 1151 11161 1117]. The horizontal lines in the left 



panel correspond to £"5 = 0, 
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the energy scales Es ^ 0, — |^7r, and others, contribute to the signal-to-noise ratio. At times 
when the nucleon correlation function is in the ground state, and the variance correlation function is 
dominated by the nucleon-antinucleon state, Es = should dominate the signal-to-noise ratio. At large 
times the variance correlation function is dominated by the 3-pion state and Es = — should 
dominate. This is modified by the finite temporal direction jTTBJ as the hadrons produced by the sources 
of the correlation function and the variance correlation function can propagate forward and backward in 
time. The calculated energy scale of the signal-to-noise ratio of the single nucleon correlation function 



obtained on 20^ x 128 anisotropic clover gauge field configurations |1151I116"1 1117j is shown in Figure 24 
It exhibits behavior that is consistent with expectations, and exceeds the long-time behavior expected 
from the Lepage argument at approximately time-slice t = 50 due to the temporal boundary conditions. 

The structure of the two-nucleon variance correlation function has significantly more structure than 
that for one nucleon. On configurations with infinite temporal extent, the proton-proton correlation 
function is of the form (neglecting interactions between the hadrons) 

(dNNit)) = 5] r^^^^(0|7V<^(x,t)7V^(y,t)iV^(O,0)iV''(O,0) |0) 

^ Z^,^e-2^-* + ..., (75) 
and the variance correlation function has the form 

= J2 Tl'^'^ri''''^^ {0\ N''{^,t)N''{y,t)N\z,t)W'{w,t) x 

x,y,z,w 

]v'(o,o)iv^(o,o)7v^(o,o)M(o,o) |o) - (eMt))' 

^ e-^^^-K (76) 
Therefore, we anticipate finding energy scales of approximately Es = 0, — ^m^ and 2Mn — Sm^ in 
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Figure 25: The energy scale of the signal-to-noise ratio, defined in Eq. (74), in the ratio of 



correlation functions that produces the shift in energy between two interacting protons and 
two isolated protons, obtained on 20^ x 128 anisotropic clover gauge field configurations |115i 
11161 I117j . The horizontal lines in the left panel correspond to Eg = 0, — and 
2Mn - 3m^. 



the signal-to-noise ratio on gauge field configurations of infinite temporal extent, see Figure 24, The 



temporal boundary conditions imposed in the present calculation introduce additional energy scales due 
to the backward propagating states. Figure [25] shows the energy scale associated with the signal-to-noise 
ratio for the ratio of correlation functions that provides the energy splitting between two interacting 
protons and two isolated protons from which the p cot 5 (p) is extracted, obtained on 20^ x 128 anisotropic 
clover gauge field configurations |115i 11161 1117j . It is clear that the energy scale of the energy splitting 
is significantly less than for the individual energies, and is consistent with zero throughout much of the 
Golden Window of time slices. This indicates that the signal-to-noise ratio associated with the energy 
splitting in the proton-proton sector, and hence the scattering parameters and bound-state energies, 
are time independent, and therefore do not degrade exponentially with time. This is an exceptionally 
important result, as it means that the extraction of NN, and more generally, multi-nucleon interactions, 
does not require an exponentially large number of calculations for each relevant correlation function. 

The signal-to-noise ratio in the SS sector is noticeably larger than in the NN and the YN sectors. 
It is likely that the improved signal-to-noise behavior in the EE sector is due to a reduced overlap of 
the source onto the multi-meson intermediate states in the variance correlation function compared to 
purely baryonic intermediate states. Such a reduction is expected based on the fact that the volume 
occupied by multiple S's is smaller than that of multiple nucleons, and serves to extend the Golden 
Window beyond its range in nucleon correlation functions. 

As discussed previously, the correlation functions for states with the quantum numbers of the triton 
(or ^He) and E^E^n have been calculated, with well-defined plateaus being observed in the EMPs in 
both channels. As the largest number of calculations have been performed in the E^E^n channel, the 
signal-to-noise ratio of its correlation function has been calculated [116j. In Figure 26 we show the 
energy-scales associated with the signal-to-noise ratios for the E^E^n correlation function as a function 
of time-slice [116j . In particular, the right-panel shows that the difference of energy-scales between 
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Figure 26: The left panel shows the energy-scales associated with the signal-to-noise ratio for 
the E^E^n correlation function [116j, as defined in Eq. (74). The horizontal line corresponds 
to rriN + 2m^ — 2m^ — ^m^^ the asymptotic energy-scale in a lattice with infinite temporal 
extent. The right panel shows the diflFerence between the signal-to-noise energy scales of 
the diagonalized EPEPti correlation function and that of the nucleon and twice that of the S 
correlation function. 



the three-baryon system and its constituents is consistent with zero within the Golden Window of 
time-slices. It is clear from this figure that the current computational technology used by NPLQCD 
would allow for the calculation of systems with A > 3 using algorithms that allow for the multi-baryon 
contractions to be performed in a reasonable amount of computer time. 



In Figure [27| the relative uncertainties in the extraction of one-, two- and three-baryon ground state 
energies obtained in Refs. |115l 11161 1117j are shown. It is somewhat remarkable from the numerical 
standpoint that the relative uncertainty per baryon is essentially independent of baryon number when 
determined within the Golden Window of time-slices. With the high-statistics calculations performed 
by the NPLQCD collaboration, an uncertainty of approximately 6 MeV per baryon has been obtained, 
which is at the threshold for performing meaningful calculations for nuclear physics (once the pion mass 
is reduced to its physical value). 

It is clear that the Golden Window of time-slices that has been uncovered in Refs. |115i 11161 1117j . 
which provides a way to evade the need for exponentially large computational resources to calculate 
multi-baryon systems, requires careful source and sink optimization. This involves two considerations 
to make optimal use of available resources: 

1. maximal overlap of the interpolating operators onto the baryon states 

2. minimum overlap of the interpolating operators onto the mesonic states in the correlation function 
dictating the variance of the baryon correlation functions 

To close this section, we make the comment that it is highly desirable to develop algorithms that 
explicitly eliminate the exponential degradation of the signal-to-noise ratio in generic correlation func- 
tions. While this is an obvious comment to make, such algorithms do not yet exist, however, there 
have been encouraging developments. It has been shown, by Liischer and Weisz [139j , that a multi-level 
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Figure 27: The relative uncertainties in the extraction of one-, two- and three-baryon ground 
state energies obtained in Refs. |115i 11161 1117j . 



scheme, involving the iterative sequential update of gauge fields and averaging the products of gauge 
links, leads to large Wilson-loops that have statistical uncertainties that are exponentially smaller than 
those resulting from the usual method of evaluation (direct evaluation of the Wilson loop on each 
member of the ensemble of gauge configurations). This multi-level scheme has been applied to the 
calculation of the Yang-Mills spectrum by inserting projectors into the transfer-matrix that restrict the 
quantum numbers of the states that can propagate forward in time |14Q[ I141j . Such projections lead to 
an exponential improvement in the time-dependence of the signal-to-noise ratio. It is conceivable that 
analogous projections can be developed for n/ = 2 + 1 calculations, for instance producing gauge field 
configurations which do not permit the propagation of single pions, thereby eliminating the long-time 
correlations in the light-quark propagators, and thereby exponentially reducing the cancellations that 
occur in forming the correlation functions for one or more baryons. There has also been progress in the 
calculation of the properties of spin-systems by using cluster-algorithms where the partition functions is 
split into the sum of contributions that individually do not suffer from the sign problem, e.g. Ref. |142j . 
However, this algorithm has yet to be transcribed into Lattice QCD calculations. 



6 Optimization through the Analysis of Simulated Calcula- 
tions 

Large computational resources are required to perform Lattice QCD calculations, especially of quantities 
that impact nuclear physics. As such, it is important to simulate the results of Lattice QCD calculations 
in order to estimate the volumes, lattice spacings and pion masses of calculations that will optimize 
the physics output with fixed resources. As an example, consider the results of a set of simulated 
calculations on an ensemble of gauge field configurations with spatial dimensions L ~ 12.3 fm at the 
physical pion mass {m^L ^ 8.7). The experimentally determined nucleon-nucleon ^S'l —^Di coupled- 
channels scattering amplitude is used to determine the two energy-levels in the finite lattice volume 
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Table 3: Simulated calculations of the shifts of the two lowest-lying energy levels of two 
nucleons in a finite lattice volume with spatial extent L ^ 12.3 fm generated from the physical 
^Si scattering amplitude (at the physical values of the light-quark masses) at varying levels 
of precision. For simplicity, the precision of the two calculated energy levels are taken to be 
equal. 





Rnnrirl Slfpifp Friprcyv i A/TpA/i 


1^^ doritirmnmTyPVpl (^A/TpA/^ 

J. V.y'WllLlllLL LLlll-LJC V I 1V±C V 1 


0% 


-3.147 


4.005 


1% 


-3.111 ±0.031 


4.015 ±0.040 


5% 


-2.95 ±0.16 


4.24 ±0.20 


10% 


-2.66 ±0.31 


3.65 ± 0.40 




Figure 28: The scattering parameters extracted from Monte-Carlo fitting to the simulated 
calculations of pcot6 shown in Table [3} The left panel corresponds to 10% precision, the 
middle panel to 5% precision and the right panel to 1% precision. 



that are below the inelastic threshold (set by |p| < using the Liischer relation in Eq. (32). In 

particular, a scattering length of o}^^^^ = 5.425 fm, an effective range of r^^^^* = 1.749 fm, and a nucleon 
mass of M^^^* = 939 MeV are used, which produce a deuteron bound by B = 2.214 MeV when d-wave 
interactions and higher order terms in the effective range expansion are ignored. Simulated results are 
then generated with 10%, 5% and 1% uncertainties by randomly distributing the centroid of the energy- 
level by the corresponding amount, and a ssig ning the corresponding uncertainty. The deuteron bound 
state energy in this finite lattice volume is E"^^ = —3.147 MeV, and the lowest lying continuum 
state appears at E^^"^ = +4.005 MeV. Given the simulated calculations shown in Table [3} the goal is to 
extract the deuteron binding energy and the scattering parameters at the varying levels of precision. 
In order to propagate the correlated uncertainties associated with the extraction of the scattering 



^^The deuteron becomes more tightly bound in the finite lattice volume due to the exclusion of momentum modes by 
the periodic boundary conditions imposed in the spatial directions. It is interesting to note that if instead of a volume 
with periodic boundary conditions, the nucleons were confined with a harmonic oscillator potential, as is common-place 
in nuclear structure calculations, the deuteron binding energy would be reduced in magnitude due to the positive energy 
contributions from the potential. 
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Figure 29: The 1-a determination of the pcotS function (shaded region) extracted from 
the simulated calculations in Table [3] (shown as points with uncertainties) . The left panel 
corresponds to 10% precision, the middle panel to 5% precision and the right panel to 1% 
precision. 




Figure 30: The deuteron binding energy extracted from Monte-Carlo fitting to the simulated 
calculations in Table [3} The left panel corresponds to 10% precision, the middle panel to 
5% precision and the right panel to 1% precision. 



parameters, we perform a Monte-Carlo analysis of this simulated set of calculations. For each value of 
and of p cot S a value is randomly drawn from a normal-distribution with mean and standard deviation 



set by the quantity and its uncertainty A two-parameter fit is performed in this simplified analysis, 
terminating the energy expansion of pcot6 at the effective range parameter, providing a pair of fit 
values to the scattering length and effective range, {a^,r^}. The distributions of these extracted pairs 



are shown in the lego plots in Figure [28j These ensembles of extractions are then used to generate ranges 
of values of pcotS for any p^ (using the functional form), from which means and standard deviations 



can be determined. The results of which are shown in Figure 29 For each pair {a^,r^}, the effective 



range expression can be used to determine the deuteron binding energy, the distribution of which is 



shown in Figure 30 , The results of the analysis of the simulated calculations are shown in Table ffl 



A number of conclusions can be drawn from the analysis of simulated calculations: 

1. Determining the shifts in the lowest two energy eigenvalues in a lattice with L ^ 12 fm with 
precision exceeding ^ 10% is sufficient to determine the deuteron binding energy with a precision 
that exceeds the 5-a-level. For the physical binding energy of the deuteron, 10% uncertainty in 
the energy shift corresponds to a ^ 0.01% uncertainty in the total energy. This is approximately 
an order of magnitude more precise than present calculations. 



It should be noted that the uncertainties associated with the simulated calculations in Table [s] are uncorrelated 
(between the and pcot 6 for each energy- level) . However, this will not be the situation in actual Lattice QCD calcula- 



tions where p cot S is determined from using Eq. ( 32 ) . For this exercise we treat them as independent (which tends to 
increase calculated uncertainties). 
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Table 4: Results of the analysis of the simulated calculations in Table [3l 





T^pii-j-pr'pv-n Ri"nrli"ncr l-i^,"nprcr"\7' i A/TpA/i 


1"1"pn "n (T T/P"ncri"n i f tti i 

kJL^ClL LCi lii^ J-JCii^Lii I ±111 J 


l~i^,fFpp1"i A/'P T-^5^"ncrp 1 f TTI 1 

J_JliCL^LlVC iLdll^C I 1111 J 


0% 


-2.214 (input) 


5.425 (input) 


1.749 (input) 


1% 


-2 180(35) 


5 447f38) 


1 737(32) 


5% 


-1.94(18) 


5.77(22) 


1.82(16) 


10% 


-1.77(35) 


5.58(51) 


1.14(42) 



2. The precision of the scattering length extraction is somewhat better than the deuteron binding 
energy, while that of the effective range is somewhat worse. This is consistent with the fact that 
the scattering length enters at LO in the momentum expansion of the scattering length, while the 
effective range enters at NLO. 

3. Having calculations on both sides of = 0, combined with the fact that the magnitude of 
the deuteron binding energy increases as the lattice volume decreases, allows for a much more 
precise determination of the scattering parameters and deuteron binding energy than would be 
the case with calculations of the deuteron binding energy alone on two lattice ensembles with 
different volumes. This is because an interpolation is required to reach the physical deuteron 
binding energy, as opposed to an extrapolation that would be required for single-level (bound- 
state) calculations in multiple lattice volumes. 

4. A NNLO analysis, including the extraction of the shape-parameter, ri requires calculating three 
energy eigenvalues in one lattice- volume, or calculations on a number of lattice volumes. Given 
the experimentally observed smallness of the shape parameters contributing to nucleon-nucleon 
scattering, high-precision calculations of the energy eigenvalues will be required for such an ex- 
traction. 

7 Other Lattice QCD Efforts in Low-Energy Nuclear Physics 

The content of this review is primarily focused on the work of the NPLQCD collaboration, and its efforts 
to extract the properties and interactions of nuclei from lattice QCD calculations. However, this agenda 
has been taken up by other collaborations also. As discussed previously, the HALQCD collaboration 
has performed both quenched and dynamical Lattice QCD calculations of baryon-baryon correlation 
functions from which non-local, energy-dependent and interpolating operator dependent baryon-baryon 
potentials are extracted [88l|90j. Unfortunately these potentials contain no more (rigorous) information 
than the location of the energy-levels in the lattice volume, and the scattering parameters that can 
be derived using the Liischer relation. Important work has been recently performed by the PACS-CS 
collaboration [133j. Exploratory quenched Lattice QCD calculations of the ^He and ^He correlation 
functions strongly suggest that both nuclei are bound for pions with a mass of ^ 800 MeV. 

An exciting development that has taken place during the last couple of years is the exploration 
of nuclei with rif = 1 Lattice QCD at strong coupling |143i I144j . In this limit there is no Yang-Mills 
contribution to the action, and instead of integrating out the fermionic fields to leave an integral over the 
gauge fields that is evaluated by Monte-Carlo, the gauge links are integrated out to leave a integration 
over the fermion fields which can also be performed in closed- form. The partition function becomes 
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a weighted sum over configurations of dimers and self- avoiding baryon loops. The nucleon-nucleon 
potential that is found, is qualitatively consistent with the modern phenomenological nucleon-nucleon 
potentials. Further, the binding energy of larger nuclei up to ^^C derived in this limit, were found to 
be qualitatively described by the first two terms in the von Weisacker formula. 

An area of Lattice QCD calculations that is peripherally related to multi-nucleon systems, but is 
having increased overlap with the study of hadronic interactions is the calculation of the excited states of 
baryons and mesons. At heavier pion masses, where many of the excited states of the nucleon and strange 
baryons are stable against strong decay, it is somewhat straightforward to calculate their spectrum with 
lattice QCD. The excited states correspond to exponentially decaying states in the correlation functions 
created by interpolating operators with the correct quantum numbers. Enormous effort has been put 
into successfully refining the basis of interpolating operators and developing techniques with which to 
cleanly separate multiple excited states from the ground-state and other excited states, for example 
Ref. [581 11451 [r46j . The situation becomes vastly more complex as the pion mass approaches the 
physical value and the resonances become unstable to strong decay. The energy eigenstates in the 
lattice volume are the baryon ground-state and the scattering states with the same quantum numbers. 
The excited states correspond to the energies for which the scattering phase-shift is at 5o = f , and the 
energy interval over which the phase-shift is "near" Sq dictates the width. Given that the location of 
the scattering states in the volume cannot be dictated in the generation of gauge field configurations, 
calculations must be performed in multiple volumes. Further, given that m^L^ 27r must be satisfied by 
calculations in order for the wavefunctions to be asymptotic at the boundary of the lattice, the location 
of multiple excited states in the lattice volume must be precisely determined, e.g. Ref. |147j . This is 
an exceptionally challenging problem, but progress is being made in extracting the properties of the 
p-meson and the A-resonance. 

An area that is not Lattice QCD, but which makes use of techniques that were developed for Lattice 
QCD calculations and impacts low-energy nuclear physics is Lattice Nuclear EFT. Low-energy nuclear 
structure and reactions are determined in the low-energy EFT using a space-time lattice, and performing 
the path integral. The conventional power-counting and perturbative expansion that are employed 
successfully for particle physics observables fail in the nuclear context as nuclei are nonperturbative 
objects. After many years of development, EFT's for nuclear physics are reasonably well developed, 
and recently progress has been made in calculating S-matrix elements in systems involving more than 
two nucleons by latticizing the EFT's ||136i I148j . 



8 Resource Requirements and the Next Decade 

A number of workshops focusing on the science need for exa-scale computing resources sponsored by 
the US Department of Energy were held during 2009. One of the workshops. Forefront Questions in 
Nuclear Science and the Role of Computing at the Extreme Scale [149j, established the need for exa-scale 
computing resources in order for the main goals of the field of nuclear physics to be accomplished. One 
of the major goals of the field that requires exa-scale computing resources is the calculation of nuclear 



forces from QCD using Lattice QCD, and Figure [31] presents an overview of current estimates of these 
requirements. 

As discussed in a previous review [97], a complete calculation of the nucleon-nucleon scattering 
amplitude, and the hyperon- nucleon and hyperon-hyperon scattering amplitudes (including multiple 
lattice spacings, volumes and light quark masses) will require sustained peta-scale resources, as shown 



in Figure |3TJ The same is true for the meson-baryon interactions. It is estimated that sustained-sub- 
peta-scale-year resources are required to perform high-precision calculations of meson-meson scattering- 
phase shifts, extrapolated to the physical pion mass (in the isospin limit), including the contributions 
from disconnected diagrams to the iso-singlet tttt channel. 
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Figure 31: Estimates of the resources required to complete calculations of importance to 
nuclear physics |149j . Except for the quantities indicated as requiring exa-scale resources, 
the resource requirements are for calculations performed in the isospin limit and without 
electroweak interactions. 



It is currently estimated that sustained peta-scale resources are required in order to calculate the 
matrix elements of electroweak operators, such as those determining neutrino-induced breakup of the 
deuteron, in the few-nucleon sector. Great progress is being made in computing single hadron ma- 
trix elements of such operators, such as the isovector axial-current matrix element in the nucleon, 
9 A jl50l I151[ I152j . While the extrapolation to the physical pion mass, and to infinite volume remain 
the subject of discussions in the community, relatively rapid progress is being made. The calculations 
of matrix elements of operators that receive contributions from disconnected diagrams remain difficult 
with currently available resources, but will be completed with peta-scale resources. A significant uncer- 
tainty in the experimentally determined properties of neutrinos comes from the uncertainties in weak 
matrix elements between nuclear states. Such uncertainties in few-nucleon systems should be reduced 
within the next decade with anticipated Lattice QCD calculations, as indicated in Figure [3T| 

Despite the first Lattice QCD calculations of three- (dynamical |116j ) and four-baryon (quenched |133j ) 
systems appearing this year, it is estimated that exa-scale computing resources will be required to ex- 
tract the nuclear interactions among three-nucleons and determine the spectrum of the a-particle. Given 
that the three-nucleon interaction is relatively imprecisely known when compared with the two-nucleon 
interactions, this calculation will have significant impact upon nuclear structure and reaction calcula- 
tions. The three-baryon interactions between strange and non-strange baryons will be calculable at the 
same time with the same level of precision with minimal additional resource requirements. 

The current discussions regarding exa-scale computing facilities suggests that it may be possible to 
see such resources deployed sometime around 2018 [149j. Clearly, such resources are required for the 
calculation of quantities of central importance to the nuclear physics program. During the next decade 
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the field will develop the ability to perform low-energy strong interaction calculations with quantifiable 
uncertainty estimation. Lattice QCD will supplement the present ability of the experimental nuclear 
physics program to determine quantities with certain precision to processes that cannot be accessed 
experimentally. 

9 Conclusion 

A central goal of the field of nuclear physics is to establish a framework with which to perform high- 
precision calculations, with quantifiable uncertainties, of strong-interaction processes occurring under a 
broad range of conditions. Quantum chromodynamics was established as the underlying theory of the 
strong interactions during the 1970's, however, nuclear physics is the regime of QCD in which its defining 
property of asymptotic freedom is hidden by the vacuum and by the phenomenon of confinement. Lat- 
tice QCD, in which the QCD path integral is evaluated numerically, is the only known way to perform 
rigorous QCD calculations of low-energy strong interaction processes. With the research and develop- 
ment into high-performance computing, nuclear physics, quantum field theory, applied mathematics, 
and numerical algorithms that has taken place over the last few decades, the field of nuclear physics 
is entering into an era in which Lattice QCD will become a quantitative tool in much the same way 
that experiments are, but with a different scope and different range of applicability. Rapid progress is 
currently being made in the calculation of the interactions among nucleons and, more generally, among 
the low-lying baryons. Present day Lattice QCD calculations are being performed at pion masses larger 
than the physical pion mass, but as exploratory calculations are now being performed at the physi- 
cal pion mass, the interactions among baryons will be known from QCD at the physical light-quark 
masses within the next several years (if computational resources devoted to these calculations continue 
to increase as they have during the last decade). 

Determining the three-body and higher-body interactions among nucleons and hyperons directly 
from QCD will be a remarkable achievement for Lattice QCD, and will provide crucial input into the 
calculations of the structure and interactions of light nuclei. During the last year the first calculations 
of three- and four- nucleon systems were reported. Such calculations are presently difficult, but are 
primarily limited by the available computational resources and not by conceptual or formal issues. It is 
anticipated that the three-nucleon interactions will be calculated with high-precision with Lattice QCD 
during the next decade. Beyond the three-body systems, we expect that some of the properties and 
interactions of light nuclei will be calculable with Lattice QCD during the same time-period. 

The dream of being able to perform reliable calculations of the interactions among multiple nucleons 
and hyperons, and of the structure and reactions of light-nuclei, directly from QCD is starting to be 
realized. The path forward is clear, and the next decade will be a truly remarkable period for nuclear 
physics. 
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